The article presents a hybrid optimization method which is implemented using the optimization toolkit DAKOTA.The method is applied in order to optimize the parameters of the recent BGC model configuration BFM17 developed by some of the papers' authors (Smith et al., 2021).



The proposed method is a combination of (many) latin hypercube samples (LHS) of the parameter space and a gradient-based optimization (gradient search, (GS)) using some best-ranked LHS samples as initial search points.

Although both, LHS and GS have been examined before in combination with population-based search methods like genetic algorithms (GAs) the proposed combination of LHS and GS seems to be new.

The application to the BFM17 BGC model is validated in a twin experiment setup where it is able to detect at least the misfit-sensitive parameters. A drawback of the twin experiment is that only a time window of 30 days is considered (according to line 294 of the paper) which does not represent the full annual cycle.

A calibration of all model parameters against 1D observations yields a new and justified parameter set for BFM17 (only a manually tuned parameter set has been provided before):

an RMSD type model-data misfit w.r.t. five tracers does significantly improve in comparison to the manually determined parameter set and, consequently, model simulation results do better agree with their observational counterparts (average mismatches drop from several standard deviations to less than one standard deviation).

However, even more convincing would be a calibration of BFM17 coupled to a global circulation model, and against global 3D observations, since (according to Smith et al., 2021) BFM17 is designated to be used in global 3D model configurations. It would also be interesting to see, how a population-based search method (alternatively incorporating LHS and GS) would perform in comparison to the proposed method. I notice that the computational demand of 25,000 function evaluations for LHS would be comparable with many iterations of, e.g., an evolutionary algorithm.

Summarizing, I suggest a minor revision to mention hybrid population-based search methods using LHS on the one hand and GS on the other hand. Also, though the main focus of the paper is on presenting and testing the proposed optimization method and providing a suitable parameter set for BFM17, I wonder if the calibration experiments reveal some general relationships between parameter values and model skills, one could elaborate on?

Some technical corrections:

- The DOI link of Vichy et al., 1998, does not work in the paper PDF.

- Line 176: "A detailed description ... has been outlined in detail ..." (2 x "detail", drop one)

- Line 305 ff: the parameter disturbance seems to be 20% instead of the stated 10%

- Line 305: The definition of S^ is difficult to read

   Perhaps replace "sensitivity metric S(pi)" by "sensitivity factor S(pi)"

   and explicitly write "S^(pi)=S(pi)/Smax with Smax=max_i(S(pi))"

- Lines 299 and 303 and Table 1: what does "nominal parameter value" mean? Is it the "baseline parameter value" from the manual calibration by Smith et al.?

- Figure 5: what is pi^ and po^? Is it the parameter values, normalized w.r.t. their credible interval length?

- Table 1: normalized RMSD values? It is according to Equation (2), i guess?

- Line 364: "year round" => "all year round"

- Line 366: "maxima" => "maximum"

The authors present a hybrid parameter estimation technique which includes a gradient-free first optimization stage and a gradient-based second stage. The two stage approach applied to a 1D coupled physical-biogeochemical ocean model appears to work well, but for a study focused on methodology, more details about it should be included.

general comments

The study makes a relatively straightforward case, the text in general and the description of the methods is mostly easy to follow, and the experiments are well motivated (maybe except for the TSE, see below). However, he manuscript frames this study as one primarily focused on the demonstration of the parameter estimation method. And here, I think, more emphasis could be placed on the parameter estimation part. (1) As a reader, it would be good to get a better idea of the computational cost in terms of runtime or number of model/function evaluations, in particular comparing the initial global search to the gradient-based second stage. (2) It seems sensible to use the two-stage (global+local) parameter estimation, but how does each stage contribute to the decrease in the cost function? Readers may ask if after 25000 simulations in stage 1, is there any need for stage 2? Or, how well does stage 2 do if started from the baseline? (3) Studies like this one often include replicate experiments, why not include one here, to show that similar parameter values are recovered when running the estimation again. (4) 25000 model evaluations are a lot, which basically prohibits the use of 3D models (BGC model coupled to 3-dimensional circulation models). What could be done to reduce runtime, would the authors consider the use of a more sophisticated gradient-free method (some are described in the introduction) instead of Latin hypercube sampling?

In the twin simulation experiment (TSE; section 4.1), it reads like all 17 state variables were used as synthetic observations. But not all 17 variables are part of the BATS or HOTS datasets, in fact, only 5 appear to be used. How do the results of the TSE change if only 5 state variables are observed? This would appear to be a much more important experiment than one using all the variables. Additionally, a little later in the section (l 312): "While these results may suggest that the least sensitive parameters could be excluded from the subsequent calibration studies, redoing the sensitivity study with our standard objective function reveals larger relative importance values for the full set of parameters.": This is a bit of an unsatisfactory result, could the discrepancy be due to including 5 compared to 17 observation types in the objective function?

specific comments

L 72: "in a high-dimensional BGC model across a range of ocean conditions, using a 1D model for vertical ocean mixing": Here it would be nice for the reader to explain a bit better what is meant by dimensionality: "high-dimensional BGC model" and "1D model for vertical ocean mixing" are used together and at this point it is unclear if "high-dimensional" refers to 3 spatial dimensions (represented at various spatial locations as 1D vertical models) or if it is referring to the dimensions of the state space. More generally, it would be useful to describe what is to follow in a bit more detail early on. Even the formulation in line 5 "simultaneous parameter estimation at multiple ocean locations" is a bit ambiguous, and the kind of model setup (1D, 3D) that is being used in the study could be mentioned earlier.

L 98: I presume \Pi is a matrix of weights and \Pi_{i,j} is a scalar weight?

L 99: "describes the misfit with observational (or other reference) data": For clarity, I suggest adding "of the model output" or similar.

Eq 2: I would have expected that the 1/\sigma term (a weight) would be contained in \Pi. In fact, the subsection makes no further reference to \Pi, and it would be useful to explain here what it is used for.

Fig. 1: It is somewhat easy to trace the line from model run and comparison with data to "Calculate model error", feeding it into DAKOTA and obtaining new parameters as output with which to run the model again. but what do the dotted gray lines mean and why does the calculated model error bypass the interface? A minor complaint is that terms like BFM17 and \partial A_j/ \partial t have not been mentioned in the text when this figure is first referenced.

L 151: "how quickly the values of J increase when the N_random simulations are sorted": That sounds a bit like J increases while the simulations are being sorted, to avoid confusion maybe use something like: "how quickly the values of J increase in the sorted N_random simulations".

L 164: "The QN algorithm reliably and efficiently converged to optimized solutions.": Is this a result from this study or a general observation? Maybe add "In our experiments ...".

L 165: "Similar to the ecosystem parameter estimation study by Matear (1995), we found that the conjugate gradient method failed to converge efficiently." It reads as if this is still meant as a justification for using QN over the conjugate gradient method. I suggest rephrasing it slightly: "In comparison, we found that the conjugate gradient method failed to converge efficiently, a result similar to that in the ecosystem parameter estimation study by Matear (1995)."

L 192: "Each CFF is represented as a vector where each element is a constituent component concentration corresponding to a state variable.": Initially, I thought that this meant that the state variables are divided up into multiple vectors according to their type/function, which is not the case according to the following paragraph. This could be explained better without adding much more text, maybe "Each CFF is a vector representing a model variable divided into the elemental constituents represented by the model, e.g., the phytoplankton CFF is a 4-element vector containing the C, N, P and chlorophyll concentration of the phytoplankton variable (see below for more details)."

L 193: "BFM17 was simplified to be a general, but computationally cheaper, model that retains the essential BGC processes for modeling a phytoplankton spring bloom ...": Is this mean in comparison to BFM56?

L 232: At this point, it is unclear from the text what parameter values are used. Interestingly, this information, plus comments about the manually adjusted parameters, is provided in the description of Fig. 4 in the next subsection, which is helpful to the reader. I would suggest moving the information that applies to both datasets (use of baseline model parameters, single-site model calibration and the multi-site calibration etc.) to the description of Fig. 3.

L 303: "nominal value": Does this refer to the optimized value?

L 304: So some parameters were only perturbed down, because perturbing them up would hit the baseline + 25% boundary? What then makes up the second perturbation case, the optimized value? This could be explained better.

L 305: "... defined as the maximum objective function evaluation between the two perturbation cases for each parameter." This sentence is difficult to understand, I suggest rephrasing it.

L 317: It would be useful to state what N_v is set to here. I presume it is 5, based on previous figures. But it is not stated directly in the text.

L 351: "improved by a factor of ... over 236 (for nitrate)": True, but this large improvement is more a function of the enormous misfit in the baseline experiment.