Optimization procedure

MODIS observations

The MODerate-resolution Imaging Spectroradiometer (MODIS) albedo product is used to optimize the albedo parameterization in ORCHIDEE. For this we use the bi-hemispherical reflectance (white-sky albedo) observation data for the years 2001-2010 at 0.5° resolution for the visible and near infrared ranges.

Additionally to that, the background albedo retrieved from the MODIS data using the Joint Research Centre Two-stream Inversion Package (JRC-TIP) is used as the bare soil albedo in ORCHIDEE instead of the soil color based scheme. These background albedo values are used as a first-guess prior estimate for \(alb_\mathsf{bg}\) and it is optimized further within the optimization procedure.

ORCHIDEE setup

The ORCHIDEE model revision 3789 is used in optimization forced by atmospheric forcing CRU-NCEP (v5.3.2 at 0.5° resolution). The reference simulation is done for 40 years (1971-2010) starting from the scratch, the last 10 years of simulation (2001-2010) are used in optimization. The ESA CCI epoch PFT map for the year 2010 (v1.6.1, with the urban class moved to PFT10, water and ice classes moved to PFT1) is used as the land cover map throughout the whole simulation.

Optimization framework

Albedo parameter optimization is done using the variational data assimilation system. It follows a Bayesian framework, in which the optimal parameter vector \(x\) can be found by minimizing the cost-function \(J(x)\), assuming that the probability distribution functions (PDFs) of the model parameter and observation uncertainties are Gaussian:

\[ J(x) = (y-H(x))^t R^{-1} (y-H(x)) + (x-x_\mathsf{prior})^t B^{-1} (x-x_\mathsf{prior}) \]

where \(y\) is the vector of observations, \(H(x)\) is the the vector of simulated outputs, \(x_\mathsf{prior}\) is the vector of prior parameter values, \(R\) and \(B\) are the prior variances/covariances matrices, describing the errors (diagonal elements) and their correlations (non-diagonal elements) in observations and parameters, respectively. The correlations are difficult to assess properly, thus they are usually neglected, and so all non-diagonal elements are set to zero. Data uncertainties (diagonal elements) are defined as the root mean squared error (RMSE) between the prior model simulations and the observations, parameters uncertainties (diagonal elements) are set to 40% of its prescribed range of variation.

The cost function is iteratively minimized, using the gradient-based L-BFGS-B algorithm (limited memory Broyden-Fletcher-Goldfarb-Shanno algorithm with bound constraints. It was developed specifically for nonlinear optimization problems and allows direct assignment of fixed bounds to each parameter. At each iteration the L-BFGS-B minimization function is called to estimate a new set of parameters. The method requires as an input the value of misfit function for the current parameter vector and its gradient:

\[ \nabla_x J(x) = 2H^t_{\mathsf{ad}} R^{-1} (H(x)-y) +2B^{-1}(x-x_b)^t \]

where \(H_\mathsf{ad}\) is the adjoint operator of the Jacobian matrix of H, having a form of \(m \times n\) matrix (\(n\) – number of optimized parameters, \(m\) – number of observation points) with each column representing a vector of the model outputs derivative versus the corresponding optimized parameter.

Using this optimization framework we have optimized two kinds of albedo parameters — \(alb_\mathsf{leaf}\) and \(alb_\mathsf{bg}\). The ranges of variation prescribed for parameters is shown in table below.

Ranges of variation for optimized parameters

PFT	2-9	10-13
\(alb_\mathsf{leaf,NIR}\)	0.15 – 0.30	0.20 – 0.40
\(alb_\mathsf{leaf,VIS}\)	0.02 – 0.10	0.04 – 0.25
\(alb_\mathsf{bg,NIR}\)	mean value ± 0.1
\(alb_\mathsf{bg,VIS}\)	mean value ± 0.1

Bare soil albedo values have been optimized individually for each grid cell, the prior value was set equal to background value from the JRCTIP package averaged through the time-series. Given the equations presented earlier, the partial derivatives of the model operator versus these parameters is expressed as:

\[ \frac{\partial H}{\partial alb_\mathsf{leaf,pft}} = \left[ 1-\exp(-\mathsf{LAI}_\mathsf{pft}) \right] \cdot frac_\mathsf{max,pft} \cdot (1 - frac_\mathsf{snow,veg}) \cdot frac_\mathsf{veg} \] \[ \frac{\partial H}{\partial alb_\mathsf{bg,ij}} = frac_\mathsf{bs,ij} \cdot (1 - frac_\mathsf{snow,veg,ij}) \cdot frac_\mathsf{veg,ij} \] where the indices \(ij\) denote the data for a corresponding pixel.

Optimization has been done in two steps explained below.

Step 1 : Leaf albedo parameters optimization

The first step primarily aims on the leaf albedo parameters optimization. To exclude the influence of the snow albedo, for this step we have filtered all the data points with snow, based on the simulation data from ORCHIDEE (if snow mass not equal to zero). Totally, it gives 403 744 (NIR range) and 403 392 (VIS range) observation data points (among 12 monthly maps at 360 x 720 grid box). Additionally to the 12 leaf albedo parameters we include in the optimization the mean background albedo (but excluding the grid cells, containing the snow throughout the whole year). Overall, we end up with 56 207 optimized parameters at this step.

Step 2 : Background albedo optimization

On the second step we fix the leaf albedo parameters to the values obtained in the first step and continue to optimize the background albedo over the full available dataset. Including back all the snow data points we have in total 709 156 (NIR range) and 708 996 (VIS range) observation data points and 61 759 optimized parameters.

Expert tuning of the snow-related parameters

The two time constants describing the snow aging and snow albedo decay have been additionally tuned during the optimization:

• the snow transformation constant (\(trans_\mathsf{snow}\)), establishing the relationship between the amount of fresh snow fall and resulting change of the snow age, was decreased from 0.3 to 0.2 m;
• the snow albedo decay rate (\(tcst_\mathsf{snow}\)), characterizing the exponential diminution of the snow albedo with its age, was increased from 5 to 10 days.

Both parameters were changed in order to counterbalance the fast aging of the snow in the ORCHIDEE model and rapid degradation of its albedo.

The albedo of the snow-covered surfaces is controlled by two PFT-dependent parameters — \(alb_\mathsf{snow,aged}\), describing the minimum albedo of the snowed surface, and \(alb_\mathsf{snow,dec}\), increasing the total albedo depending on the snow age, so that the maximum albedo equals to \(alb_\mathsf{snow,aged} + alb_\mathsf{snow,dec}\). These parameters, being obviously too low in the default ORCHIDEE configuration, have been also manually tuned in order to better match the MODIS albedo data during the snowing time. This kind of tuning was done by controlling the amplitude and dynamics of the simulated albedo on the regional seasonal cycle graphs (shown in the Results section).