2.4.1 The principle of minimisation

The optimisation methodology used in this study relies on a Bayesian statistical formalism (see for instance Tarantola, 2005, Chapter 3). Assuming that the errors associated with the parameters, the observations and the model outputs follow Gaussian distributions, the optimal parameter set corresponds to the minimum of a cost function, J(x), that measures the mismatch between (i) the observations (y) and the corresponding model outputs, H(x), (where H is the model operator), and (ii) the a priori (x_b) and optimised parameters (x), weighted by their error covariance matrices:

R and P_b are the prior covariance matrices for the observation and parameter errors, respectively. The error covariance matrices are difficult to assess and thus neglected here, as in most studies (i.e. R and P_b are diagonal). The observation errors (variances) have been defined as the mean-squared difference between the observations and the prior model simulations. This reflects not only the measurement errors but also possibly includes significant model errors. The exact error values assumed for NEE and LE for each site can be found in Tables S3–S4.

As explained in the introduction, two classes of method to minimise the function J(x) are compared; both are described below. Note that using non-Gaussian errors would significantly complicate the application of one class (gradient-based) and is thus out of the scope of this study. MacBean et al. (2016) examined the issues that may arise when using Gaussian assumptions in gradient-based minimisation algorithms; however, they found that the algorithm used in this study could account for quasi non-linearity. Moreover, in the case of ORCHIDEE we have shown in an earlier study (Santaren et al., 2007) that most parameter errors follow Gaussian distributions.

2.4.2 Gradient-based minimisation algorithm

For the gradient-based approach we chose the quasi-Newton algorithm L-BFGS-B to iteratively minimise the cost function (limited memory Broyden–Fletcher–Goldfarb–Shanno algorithm with bound constraints; see Byrd et al., 1995), referred to as BFGS. It was specifically developed for quasi-non-linear optimisation problems and allows for direct assignment of bounds to each parameter.

Each iteration requires the evaluation of the cost function and its gradient with respect to each parameter. Generally, the gradient can be calculated with a finite-difference approximation (i.e. the ratio of the change in model output to the change in model parameter). However, it can be done more precisely using the tangent linear (TL) model (linear derivative of the forward model). For ORCHIDEE the corresponding TL model has been derived with the automatic differentiation tool TAF (Transformation of Algorithms in Fortran; see Giering et al., 2005), following code cleaning and structural adjustments (without changing the physics) to allow the use of TAF (a significant effort that took around 2 years). It was used to calculate the gradient for all parameters except for two phenology parameters, K_pheno,crit and C_Tsen, associated with threshold functions. For these two parameters, the finite-difference method is more appropriate and was used here.

The search is terminated when the relative change in the cost function becomes smaller than 10⁻⁴ during five consecutive iterations. On average the iterative process stopped within 30 iterations. Given that the computation of a TL model run is about twice as long as the standard run of ORCHIDEE, the total computation time for one BFGS optimisation for one site with 30 parameters is equivalent to around 1800 standard model runs.

Gradient-based algorithms have a potential drawback with non-linear models as they are more likely to converge to a local minima of the cost function instead of the global one. In order to assess the extent of this problem we performed a set of independent assimilations starting with 16 different random initial first guesses for the parameter vector, as in Santaren et al. (2014) (see Sect. 2.5).

2.4.3 Random search minimisation algorithm

For the random search approach, we chose the genetic algorithm (GA). The GA is a meta-heuristic optimisation algorithm, belonging to a larger class of evolutionary algorithms that follows the principles of genetics and natural selection (Goldberg, 1989; Haupt and Haupt, 2004). It considers the vector of parameters as a chromosome with each gene corresponding to a given parameter. The algorithm works iteratively, filling a pool of n chromosomes at each iteration. The starting pool is created with randomly perturbed parameters. The chromosomes at each following iteration are created from the randomly chosen chromosomes of the previous iteration by one of the two processes:

• exchange of the gene sequences of two parental chromosomes (“crossover” process),

• random perturbation of the selected genes of one parental chromosome (“mutation” process).

The resulting pool is then filled with n best chromosomes from both parental and offspring pools, corresponding to the n lowest cost function values. Furthermore, at each step the chromosomes are ranked in ascending order of the corresponding cost function values and random selection of the new parents is weighted by those ranks in order to guarantee that the best chromosomes produce more descendants (“selection” principle).

The GA performance is sensitive to its particular configuration. In this study, we use the same GA configuration as in Santaren et al. (2014) who tested different GA configurations and chose one, based on computational efficiency (to locate the cost function minimum):

• number of chromosomes in the pool is 30;

• number of iterations is 40;

• crossover/mutation ratio is 4 : 1;

• number of gene blocks exchanged during crossover is 2;

• number of genes perturbed during mutation is 1.

For a better convergence of the cost function we reduce the ranges of parameter variations after the 30th iteration by a factor of 0.25 with the centre point at the current best-guess value.

The GA does not require any gradient calculations; therefore, one chromosome rank estimate requires one standard model run and the total computational time for the GA assimilation for one site equals 1200 runs. Contrary to gradient-based methods, random search algorithms allow us to use any form of PDFs for the observation and parameter uncertainties, and thus to use non-Gaussian PDFs.

2.4.4 Posterior uncertainty

Assuming Gaussian prior errors and linearity of the model in the vicinity of the solution, the posterior error covariance matrix of the parameters, A, can be approximated by

where H is the model sensitivity (Jacobian) at the minimum of J(x) (see for instance Tarantola, 1987). From A we can compute error correlations between parameters, a diagnostic that is used to evaluate the information content brought by the observations to discriminate each parameter.

3.1.1 Single site optimisation: comparative performance of the methods (S_BFGS vs. S_Genetic)

The first objective of the study is to assess which method performs best for a SS optimisation. Figure 2a and b compares the performance of the two methods for all sites in terms of model–data RMSD reduction for the NEE fluxes. The similar comparison for the LE fluxes can be found in the Supplement (Fig. S2). The RMSD reductions are expressed in percentage; they correspond to the median value (left plots) and the 5th–95th percentile range (R₉₀; right plots) across the 16 random first-guess tests.

At all sites of all PFTs the SS GA (S_Genetic) provides better fit to the data (median across all 16 first-guess tests) than the SS gradient-based algorithm (S_BFGS), with only a few exceptions (see Figs. S1 and 2a). The S_Genetic leads up to 40 % of additional model–data RMSD reduction compared to the S_BFGS algorithm, with a mean improvement across all sites of 10 %. In most cases the worst posterior RMSD value achieved with the S_Genetic algorithm is better than the median value for the S_BFGS method. The S_BFGS approach is slightly better only at eight sites with no more than 4 % of additional RMSD reduction.

The spread of the results obtained with the use of different first-guess parameters (16 tests) is also much lower for the S_Genetic compared to the S_BFGS method, for nearly all sites (Fig. 2b). On average the R₉₀ values are 3.6 times lower with the GA – the S_Genetic method is only slightly disadvantageous for 3 sites out of 78, although the spread across the 16 first guesses remains low (<10 %). Overall, the mean spread of the RMSD reduction for the S_Genetic case is around 10 %. This clearly indicates that the GA, following the set up described in Sect. 2.4.3, is able to find nearly the same minimum of the cost function independently from the choice of the first-guess parameters (with a significant improvement of the model–data misfit, see Fig. 1). On the contrary, the results of the S_BFGS method are highly dependent on the first-guess model parameters, with a spread above 22 % for half of the sites. Finally, if we compare the best achieved RMSD reduction within the 16 first-guess tests (see Tables S3 and S4 in the Supplement), we obtain similar performances for most of the sites of 5 PFTs out of 7, excluding TropEBF and BorDBF, where the S_Genetic method still produces better results than the S_BFGS method.

Figure 2Comparison of the performances between the gradient-based and genetic methods in terms of model–data RMSD reduction (%) obtained for the daily NEE fluxes. The left panels show the medians from the results obtained within 16 optimisation tests with random first-guess parameter values for each site; the right panels show the spreads between the 5th and 95th percentiles (R₉₀) of the same distributions.

Download

With the gradient-based method (S_BFGS), using the standard ORCHIDEE model parameters as a first guess does not guarantee an optimal reduction of the cost function – the corresponding posterior RMSD could be either the lowest or the highest one of the 16 tests (these cases are shown by a circle, as opposed to a cross, in Fig. S1 that displays the posterior RMSD for all sites and all 16 tests). Although we have used the same random parameter sets for each site of a given PFT, the “optimal” first-guess parameter set (i.e. the one providing the largest cost function reduction) differs between sites. It indicates that the shape of the cost function varies between sites and that there is no optimal prior first guess with the S_BFGS gradient method, which is prone to getting caught in local-minima if the assumption of model quasi-linearity is not met. Overall, this highlights one weakness of gradient-based methods and the need to perform several independent assimilations starting from different first-guess parameter values (see Sect. 3.1.4).

3.1.2 Multi-site optimisation: comparative performances of the methods (M_BFGS vs. M_Genetic)

The reduction in NEE RMSD for the MS optimisations (M_BFGS for the gradient-based method and M_Genetic for the Genetic algorithm) is illustrated in Fig. 2c, d. First, the MS optimisations provide lower RMSD reduction than the SS optimisations (Fig. 2a, b vs. c, d). This is the trade-off between fitting a specific site versus fitting an ensemble of sites representing more accurately the diversity of plant, soil and climate for a given PFT. We then notice that for few sites the RMSD is increasing after the optimisation (i.e. negative value of the RMSD reduction): 1 site for TropEBF and C₃ grass PFTs and 3 sites for TempENF and BorENF. On average these sites have only 1 or 2 years of data with large prior observation errors and thus a smaller weight in the cost function compared to the other sites of the same PFT. This could thus explain that the optimisation degrades the model-data fit at these sites but it could also indicate that there is no optimal parameter set improving the model-data fit at all sites, suggesting the need to refine the PFT description and/or to improve the model structure (see a list of potential impacts of missing processes in the ORCHIDEE model in Sect. 2.1). Note that MS optimisations are likely to be more impacted by model structural errors than SS optimisation, as these errors may have different impacts at each site. Limits arise when the objective is to study the response of a specific ecosystem (and not a generic PFT) to external drivers as the MS parameter set might be suboptimal.

The performance of the two methods differs between the PFTs. For TropEBF, TempENF and C₃ grass PFTs the GA still provides on average significantly larger RMSD reduction than the gradient method, with ratios between the two of 2.0, 1.8 and 1.4, respectively. For the other four PFTs both methods show much smaller differences and for TempDBF the average RMSD reduction is exactly the same. Overall, at 50 sites (from 78) the M_BFGS method provides performances comparable to the M_Genetic method (RMSD reductions differ by less than 25 %) or even slightly better (see Tables S3/S4 in the Supplement for all numbers) and the mean additional RMSD reduction for the M_Genetic method is only 6 %.

The spread across random first guesses between the two methods is much more comparable than for the SS optimisations. The R₉₀ values for the gradient-based method are only considerably larger for half of the sites than the GA method (with values up to 2 times larger), except for two PFTs (TempENF and TempDBF) where they are similar on average (Fig. 2d). This suggests that increasing the number of observations and/or capturing a greater range of sensitivity to the parameters acts to linearise or smooth the cost function, thereby ensuring that the gradient-based method does not get as caught in local minima as in the SS optimisations.

The last point to notice is that using the standard ORCHIDEE parameter set as a first guess with the M_BFGS method always leads to significant improvement of the model-data fit (see Fig. S1) with RMSD reduction at the same level as the median RMSD reduction for the M_Genetic method. This was not the case for the SS optimisation (see Sect. 3.1.1). With a MS optimisation the gradient-based scheme is less dependent on the first guess, likely due to the smoothing of the cost function as discussed above and therefore performing several random first-guess tests is less needed. Overall for the MS optimisations the choice of the minimisation algorithm seems thus less crucial than for the SS cases.

Nonetheless, the MS method conceals certain limitations. The optimisation may not perform well (i.e. not lead to a large minimisation of the cost function) if the different sites have very different behaviours in terms of the carbon and water cycle responses to climate forcing. This case points to the need to reconsider the PFT geographical description with possible further regional split; it is suggested for TropEBF and C₃ grass PFTs. Additionally, MS optimisations, that requires large computing time, are more complicated to set up, with the need to have coherent observation errors between sites. We need to avoid that one or few sites dominate the cost function because of a too low observation error (measurement and model errors grouped in the R term) and thus prevent the optimisation to fit all the other sites.

Figure 3Model–data RMSD reduction (in %) obtained for the NEE fluxes as a function of the number of runs performed with random first-guess parameter values (for each configuration, S_BFGS, S_Genetic, M_BFGS and M_Genetic). For each number of first guesses (x axis) all possible combinations across the 16 optimisation tests (i.e. 16 first guesses) are considered and the maximum RMSD reductions are calculated; the mean of these maxima is reported on the y axis. The results are averaged across all sites for each configuration.

Download

3.1.3 Benefit of multiple first guesses for a gradient-based approach

We now investigate more quantitatively the benefit of using several first-guess tests (16), especially for the gradient-based algorithm (see Sect. 2.5).

In order to analyse the influence of the number of first-guess tests on the optimisation performance, for each configuration we bin the 16 first-guess optimisations in all possible groups of n elements and defined for each group the maximum RMSD reduction. We then take the mean values of all these maxima. The results are shown in Fig. 3 for all configurations (S_BFGS, S_Genetic, M_BFGS and M_Genetic) as a function of the number of first guesses (n) for the NEE fluxes (with the similar results in Appendix, Fig. S3 for the LE fluxes). Only the mean across all PFTs is displayed and analysed. This illustrates the gain of adding more first-guess tests in the different configurations.

On average for the SS configurations, the GA (S_Genetic) with any number of first guesses leads to 50.5–55 % of RMSD reduction for the NEE fluxes; the increase in RMSD reduction between 1 and 16 first guesses is small (≈4 %). However, the gradient-based method (S_BFGS) needs 15 independent first-guess optimisations in order to reach the same level of RMSD reduction as with one first guess of S_Genetic and there is a substantial increase in the RMSD reduction from one to five first guesses (≈9 %), while the improvement afterward is much smaller. For the S_Genetic method, the increase in the RMSD reduction is only significant from one to two first-guess tests.

For the MS configurations, the performances of the gradient-based method are more similar to that of the GA (although both with smaller RMSD reduction, as discussed in the previous section). On average, the M_BFGS method needs three first-guess tests to reach the performance of the M_Genetic method obtained with one test. Like with the SS configuration the increase in the RMSD reduction is substantial for the M_BFGS method between 1 and 5 first guesses (≈10 %) and much smaller afterwards.

Overall this analysis highlights that using several first guesses is crucial for the gradient-based method and that at least 3 to 5 tests are required to obtain a RMSD within 5 % of the optimal value obtained with 16 tests (for both S_BFGS and M_BFGS). The GA is much less sensitive to the number of first-guess tests performed. Using a MS configuration reduces the difference between the two algorithms and requires less random first-guess tests for the M_BFGS method to achieve performances comparable to the M_Genetic method. Note finally that increasing the number of first-guess tests to much larger values would lead to a convergence of the method performances.

3.2.1 Estimating the correct posterior parameter value – pseudo-observation experiments

Figure 4 shows the estimated parameter values, the standard prior values and the “true” values for the pseudo-observation experiments with TempDBF. We display the mean and standard deviation across the 16 first-guess tests and for the SS assimilations we take the mean across all sites.

On average both BFGS and GA are not able to retrieve exactly the true value for all parameters (i.e. values used to generate the pseudo data), although many of the parameters are relatively well estimated by the optimisation schemes. Out of 28 parameters, on average across all methods and all first guesses, 14 have posterior values that are within 5 % of the true value, using the physical parameter range to define the percentage (i.e. true value ±5 % of the allowed range of variation); 6 parameters are between 5 and 10 %; 7 between 10 and 20 %, and 1 parameter (Z_crit,litter) is 29 % from the true value. These differences are mainly related to differences in the sensitivity of the model outputs (NEE and LE) with respect to the different parameters (with lower sensitivity for Z_crit,litter). The most constrained parameters (5 % range with the small deviation around true value) are G_s,slope, c_Topt, K_pheno,crit, C_Tsen, L_age,crit, τ_leafinit, K_soilC and Q₁₀. We can thus expect more robust results with these parameters using real data. However, the uncertainties associated to the real data, and the fact that these uncertainties may not be adequately described in the observation error covariance matrix, will complicate the parameter retrieval (see MacBean et al., 2016). Note that for many other parameters the error bars (standard deviation across all first-guess tests) encompass the true values (V_cmax, c_Tmin, c_Tmax, SLA, K_GR, MR_offset, MR_slope, HR_Hc, HR_Hmin, Z_decomp, K_albedo,veg, K_rsoil).

As it was already outlined in Santaren et al. (2014), the differences obtained between the true and optimised parameter values are likely due to equifinality problems (i.e. multiple solutions achieving the same overall global minimum of the cost function) or the existence of local minima where the algorithm is trapped. These issues vary between parameters and are related to the sensitivity of the model outputs to each parameter, the prior parameter uncertainty and the level of error correlation between parameters. Existing correlations and anti-correlations between the impacts of different parameters prevents from retrieving all of them with the chosen set of observations (daily means of NEE and LE) and the optimisation frameworks that are used.

Differences between the methods exist but are not systematic. There are a few parameters that perform well only with one of the methods, both for SS and MS optimisations. For example, c_Topt, τ_leafinit and Z_decomp are correctly estimated with only the gradient-based approach. On the contrary, the generic algorithm performs much better for parameters like LAI_max, SLA, L_age,crit, K_LAIhappy and K_rsoil. A few parameters are not well estimated by both methods, like F_stressh, HR_Hb and Z_crit,litter, having the largest difference with respect to the true value. The estimated errors (Eq. 2) associated to these poorly retrieved parameters are relatively large, and for F_stressh and HR_Hb they are highly correlated with other parameters; we should thus be very cautious when interpreting their value with real data optimisations. Overall, we observe that 19 parameters out of 28 are better estimated on average by the GA (S_Genetic, M_Genetic) than by the gradient-based method (S_BFGS, M_BFGS). This is coherent with the results of the fit to the data (Sect. 3.1), indicating that the gradient method is more likely stuck in local minima than the GA.

Figure 5Mean posterior model parameter values for temperate deciduous broadleaf forest PFT. Results for the four different methods (S_BFGS, SB; S_Genetic, SG; M_BFGS, MB; and M_Genetic, MG) are displayed: mean value and standard deviation across 16 random first-guess tests. Grey horizontal lines represent prior values, which are equal to the default ORCHIDEE model values. The vertical axis limits represent the range of parameter variation.

Download

3.2.2 Spread in parameter values across methods and first-guess tests – real data experiments

Using real data, Fig. 5 displays the mean posterior estimates of the TempDBF parameters across all first-guess tests for the four optimisation methods (S_BFGS, S_Genetic, M_BFGS and M_Genetic). For most of the parameters the mean optimised values are relatively similar between the BFGS and the GAs, with some exceptions. Twelve parameters (G_s,slope, c_Topt, c_Tmin, K_pheno,crit, C_Tsen, LAI_max, K_LAIhappy, Q₁₀, MR_offset, K_GR, HR_Hc, K_z0) out of 28 show posterior differences between S_Genetic and S_BFGS methods that are lower than 5 % of the physical range for each parameter, while 7 (c_Tmax, F_stressh, K_wroot, SLA, τ_leafinit, HR_Hmin, K_albedo,veg) show differences between 10 and 20 %, whereas K_rsoil goes up to 36 %. For M_Genetic and M_BFGS methods, only 7 parameters are within the 5 % variation range (G_s,slope, C_Tsen, LAI_max, K_z0 as for the SS approaches, plus HR_Ha, Z_decomp, Z_crit,litter), and 14 parameters have differences over 10 % (same list as for the SS methods excluding c_Tmax, plus c_Topt, c_Tmin, K_pheno,crit, L_age,crit, K_LAIhappy, K_soilC and MR_slope). Note that most parameters that are well constrained (i.e. small spread across the optimisation methods) also correspond to the most constraint ones in the pseudo-observation experiments (like G_s,slope, c_Topt, K_pheno,crit, C_Tsen, Q₁₀, HR_Hc). However, few parameters demonstrate different behaviour between the real and pseudo-data tests (like K_LAIhappy, MR_offset, LAI_max).

For certain parameters (c_Topt, LAI_max, SLA, L_age,crit, K_LAIhappy, K_GR, K_rsoil) we obtain significant differences in the estimated values between BFGS and GA, for both SS and MS cases. Whereas the gradient-based methods looks for the optimal parameter set in the vicinity of the first-guess set-up, the random search algorithm may jump to a completely different parameter state during one iteration. Given that with pseudo-data the GA manages to find the true values for these parameters much more precisely than the BFGS, we can speculate that the GA provides more optimal posterior estimates in the real data experiments as well.

Another feature is that for the SS and the MS cases, both algorithms (BFGS and GA) lead, for most parameters, to a similar dispersion of the posterior estimates from the ensemble of 16 first guesses (see Fig. 5). However, for some parameters (F_stressh, SLA, MR_offset, MR_slope, HR_Ha, HR_Hb, HR_Hc, HR_Hmin, Z_crit,litter, K_albedo,veg, K_rsoil) the GA method gives, surprisingly, much higher distribution of the posterior values. It corresponds to the parameters that have also not been correctly estimated in the pseudo-observation tests (see Sect. 3.2.1) with estimated value outside the 10 % range around the true value. Although not intuitive, the random sampling over the parameter space with the GA could explain that for each first guess the GA method may end up exploring a larger domain of the parameter space and thus converge to more different parameter sets than the gradient-based method.

For each parameter, the dispersion obtained with the different first guesses is lower for the MS case compared to the SS case, for both algorithms. On average the mean dispersion for the SS approaches is 25 % of the parameter range (with minimum and maximum spreads of 7.3 % and 42 %, respectively), and for the MS approaches it is only 14 % (between 3.4 % and 32 %). This reflects that with the MS approach the shape of the cost function is more smooth due to the combination of different NEE/LE responses to the parameter variations (differences induce mainly by different climate, soil type and plant species), which in turns facilitate the convergence to the global minimum. For the SS case, the shape of the cost function may significantly differ between sites leading to a larger parameter spread in response to multiple first guesses (compared to the MS case). This is specifically illustrated with the parameters K_wroot, C_Tsen, Q₁₀, HR_Ha, HR_Hb and K_z0. The only exception is for K_soilC, a parameter that is specific to each site even in the MS approach (see Sect. 2.2).