Parameter calibration and stomatal conductance formulation comparison for boreal forests with adaptive population importance sampler in the land surface model JSBACH

We calibrated the JSBACH model with six different stomatal conductance formulations using measurements from 10 FLUXNET coniferous evergreen sites in the boreal zone. The parameter posterior distributions were generated by the adaptive population importance sampler (APIS); then the optimal values were estimated by a simple stochastic optimisation algorithm. The model was constrained with in situ observations of evapotranspiration (ET) and gross primary production (GPP). We identified the key parameters in the calibration process. These parameters control the soil moisture stress function and the overall rate of carbon fixation.

The JSBACH model was also modified to use a delayed effect of temperature for photosynthetic activity in spring. This modification enabled the model to correctly reproduce the springtime increase in GPP for all conifer sites used in this study. Overall, the calibration and model modifications improved the coefficient of determination and the model bias for GPP with all stomatal conductance formulations. However, only the coefficient of determination was clearly improved for ET. The optimisation resulted in best performance by the Bethy, Ball–Berry, and the Friend and Kiang stomatal conductance models.

We also optimised the model during a drought event at a Finnish Scots pine forest site. This optimisation improved the model behaviour but resulted in significant changes to the parameter values except for the unified stomatal optimisation model (USO). Interestingly, the USO demonstrated the best performance during this event.

1 Introduction

Plants exchange carbon dioxide (CO₂) and water vapour (H₂O) with the atmosphere. Sufficient soil water, irradiance and adequate temperature are required to maintain the exchange rates during the growing season. Disturbances in these conditions such as drought, cold temperature or low radiation cause the plants to respond to the environmental stress via stomatal closure and the decrease in photosynthesis and transpiration (Lagergren and Lindroth, 2002; Mäkelä et al., 2004; Gao et al., 2017). The capability of plants to recover from such events depends on species and their adaptation to site conditions (Kozlowski and Pallardy, 2002). Stress is part of the normal annual cycle of the plants, but occasionally it may exceed the limits of recovery.

Soil water deficit and high water vapour pressure deficit can result in suppressed plant transpiration (Bréda et al., 1993; Kropp et al., 2017). Globally, soil drought has been recognised as one of the main limiting factors for plant photosynthesis (Nemani et al., 2003), and boreal forests are known to occasionally suffer from soil drought (Muukkonen et al., 2015; Gao et al., 2016). The recovery of photosynthetic capacity in spring has been connected to temperature history and to the frequency of severe night frosts (Bergh et al., 1998; Bergh and Linder, 1999), which can reverse the recovery. Understanding and correctly modelling these phenomena are especially important for boreal forests (Bonan, 2008) under changing environmental conditions.

Ecosystem and land surface models, describing the plant photosynthesis, transpiration and soil-hydrology-related processes, usually include descriptions and parameterisations for various stress effects. These parameters often lack a theoretical foundation (Gao et al., 2002; Medlyn et al., 2011), and descriptions of vegetation drought response and phenology have been recognised to need better formulations and design (Richardson et al., 2012; Powell et al., 2013; Xu et al., 2013; Medlyn et al., 2016). These deficiencies restrict a model’s predictive capability under changing environmental conditions and call for specific parameterisations for different plant types and vegetation zones.

Stomatal conductance models describe the pathway of CO₂ and water through the leaf stomata by an electric circuit analogy (Nobel, 1999). The variations in stomatal opening and mesophyll structure are interpreted as resistances to water flow and the process is idealised via generalised parameterisation. Stomatal conductance models mainly differ in their choice of variables driving the stomatal closure, and their performance has recently been assessed in modelling studies by, e.g., Egea et al. (2011), Knauer et al. (2015) and Franks et al. (2018). However, it can be hypothesised that the choice of the stomatal conductance model affects the ecosystem model parameters more broadly as the stomatal conductance formulations vary in their responses to the different conditions. A holistic assessment of the performance of the stomatal conductance models together with ecosystem model parameter optimisation has been missing.

In many other studies, where the aim has been to optimise land surface model parameters, the optimisation is based on estimating the gradient of the cost function: Knorr et al. (2010) for JSBACH, Kuppel et al. (2012) and Peylin et al. (2016) for ORCHIDEE, and Raoult et al. (2016) for JULES. Gradient-based methods are faster than Markov chain Monte Carlo (MCMC) methods as they strongly steer the sampling process to reach a minimum in the cost function (see, e.g., Gelman et al., 2013). This approach also enables a more indefinite setting of parameter ranges (limits for acceptable parameter values) when compared to methods that sample the full parameter space. However, they are prone to get stuck in local minima, especially when the dimensionality of the parameter space increases. In the last few years, similar parameter estimations have also been done for CLM (Community Land Model) by Post et al. (2017) using the DREAM_(zs) (MCMC) algorithm with multiple chains and for JULES by Iwema et al. (2017) with the BORG algorithm that employs multiple optimisation algorithms simultaneously. The DREAM algorithm is fully iterative, which limits the number of parallel processes to the number of parallel chains in use (when we do not account for the possibility of the model parallelisation, that can be substantial). The applicability of the BORG algorithm is dependent on the algorithms in use and the expertise of the user (to choose the right algorithms, etc.).

APIS (adaptive population importance sampler) is a Monte Carlo (MC) method that can be run iteratively as presented by Martino et al. (2015), but it is also straightforward to parallelise, since all samples prior to each adaptation (in our simulation 2000 draws) can be drawn and estimated simultaneously. This latter feature is useful to decrease the amount of real time required to run the algorithm when computer resources are not the limiting factor – APIS requires considerably fewer sequential estimates than typical Markov chain methods. In the iterative mode, automatic stopping rules can be easily implemented to indicate when additional samples are not required to improve the estimates. The APIS algorithm samples the full parameter space (as do MCMC methods) and can utilise a mixture of parameter prior distributions. Therefore, APIS can estimate complicated multidimensional probability distributions with relative ease. These aspects make APIS an attractive alternative to the other sampling and optimisation methods mentioned above.

In this study we apply the land surface model JSBACH for 10 boreal coniferous evergreen forest eddy covariance sites to examine the performance of different stomatal conductance models, and their effect on calibrated parameters related to photosynthesis, phenology and hydrology. First, we utilise APIS to sample the full parameter space with the different stomatal conductance formulations and to locate different modes of the target distributions (peaks of high probability). Second, using the distributions generated by APIS as the prior distributions, we optimise the parameters using a simple stochastic optimisation method. Finally, we assess the inter-site variability and the robustness of the calibrated parameters together with different stomatal conductance formulations. Optimised parameters for a specific drought are also investigated and compared with the parameters for the general optimisation.

2 Materials and methods

We will next introduce the measurement sites, followed by the model and modifications made to it. Afterwards we will give a general overview of the simulations as well as the sampling process, the algorithms and methods used to analyse the results.

2.1 Sites and measurements

We use data from 10 FLUXNET (Baldocchi et al., 2011) sites characterised as coniferous evergreen forests. Site descriptions with appropriate references are provided in Table 1. The site-level half-hourly eddy covariance (EC) measurements were quality checked and gap-filled when needed to produce continuous half-hourly and daily time series. The gap-filled and low-quality (based on FLUXNET data quality flags) measurements were masked, and the daily aggregates (usually means) were accepted as part of the calibration process if at least 60 % of the values between 04:00 and 20:00 (i.e. daytime measurements) for that day were unmasked. The daily aggregates of evapotranspiration (ET) and gross primary production (GPP) were used to calibrate and validate the model, whereas the half-hourly data were used as climate forcing (as explained later in Sect. 2.4.

Based on the quality and quantity of their respective measurements, the sites were divided into calibration and validation sites. Essentially, if we have enough data from a site, it is used for both calibration and validation purposes. We required the site to have at least 8 years of measurements, where the first five were used for calibration, and the consecutive three for validation. Otherwise we used the site only for a 3-year validation. The FLUXNET datasets were missing both the long- and shortwave radiation for the two Russian sites – Fyodorovskoye (RU-Fyo) and Zotino (RU-Zot). These were generated from ERA Interim data. The soil types of all of these sites can mostly be identified as mineral soils with varying sand, clay and peat contents. Fyodorovskoye and Poker Flat (US-Prr) are natural peatlands, and Lettosuo (FI-Let) is a drained peatland site.

The measurement error in the EC flux data was separated into systematic and random errors. The main systematic errors (density fluctuations, high-frequency losses, calibration issues) were taken into account as part of the post-processing of the data, and the random errors tend to dominate the uncertainty of the instantaneous fluxes. The random error is often assumed to be Gaussian but can be more accurately approximated by a symmetric exponential distribution (Richardson et al., 2006). It increases linearly with the magnitude of the flux, with a standard deviation typically less than 20 % of the flux (Richardson et al., 2008; Rannik et al., 2016). Our treatment of the measurement (and model) errors is explained in Sect. 2.9.

Chen et al. (2006)Chen et al. (2006)Kolari et al. (2009)Aurela et al. (2015)Thum et al. (2007)Launiainen et al. (2016)Chen et al. (2006)Launiainen et al. (2016)Kelliher et al. (1998)Ikawa et al. (2015)

Table 1 Descriptions for the sites used in this study sorted by their FLUXNET identifier. The first six sites are used for both calibration and validation purposes, with the first 5 years of each site used for calibration. The last 3 years as well as the last four sites are used for validation only. The reported elevation is in metres above sea level, LAI is the one-sided leaf area index, and the average stand age is in years, along with average annual precipitation (P) in millimetres and temperature (T) in degrees Celsius.

2.2 The JSBACH model

JSBACH (Kaminski et al., 2013) is a process-based ecosystem model and the land surface component of the Earth System model of the Max Planck Institute for Meteorology (MPI-ESM). We ran JSBACH offline using meteorological measurements from the flux towers to force the model. Implications of this one-way coupling with the atmosphere include lack of feedback from the surface energy balance to the atmosphere; i.e. latent and sensible heat fluxes and surface thermal radiation do not directly affect prescribed air temperature or humidity. Similarly, the feedback of the surface to the vertical transfer coefficients within the atmospheric surface layer is missing as the wind speed that drives mixing is prescribed. Furthermore, since we use site-level data (each site is represented as a single grid point), the grid resolution does not affect the results.

We focus only on the most essential parts of JSBACH relating to our work. A more complete model description with details on, e.g., soil heat transfer, water balance and coupling to the atmosphere can be found in Roeckner et al. (2003), whereas Raddatz et al. (2007) provides a more descriptive synopsis on land–surface interactions, Reick et al. (2013) complements both with an addition of land cover change processes, and Hagemann and Stacke (2015) introduces soil hydrological mechanisms within a multilayer scheme applying five layers.

In JSBACH, the land surface is divided into grid cells, which are split into bare soil and vegetative areas. The vegetative area is further divided into tiles representing the most prevalent vegetation classes, called plant functional types (PFTs) (Reick et al., 2013). In our site-level simulations, the model was set to use only one PFT: coniferous evergreen trees. The seasonal development of leaf area index (LAI) for the trees is regulated by air temperature and soil moisture with a single limiting value (for all sites) for the maximum of LAI. This maximum value was fixed and the site-specific fractions of vegetative area were adjusted to reproduce the measured site-level LAI.

The predictions of phenology are produced by the Logistic Growth Phenology (LoGro-P) sub-model in JSBACH (Böttcher et al., 2016). Photosynthesis is described by the biochemical photosynthesis model (Farquhar et al., 1980). Following Kattge et al. (2009), we set the maximum electron transport rate (J_max) at 25 ^∘C to 1.9 times the maximum carboxylation rate (V_C,max), which is in line with, e.g., Leuning (2002) and Ueyama et al. (2016). The photosynthetic rate is dependent on the stomatal conductance formulation used, introduced in Sect. 2.3. Radiation absorption is estimated by a two stream approximation within a three-layer canopy (Sellers, 1985). Especially in sparse canopies, radiation absorption is affected by clumping of the leaves which is here taken into account according to the formulation by Knorr (1997).

Parameters detailing site-specific soil properties, such as soil porosity and field capacity, were derived from FLUXNET datasets and the references in Table 1. We approximated the soil composition and generated these properties following Hagemann and Stacke (2015).

2.3 Modifications to the JSBACH model

All parameters of interest, presented in Table 2, were extracted from the JSBACH model code to an external file to facilitate the simulations. The default values of newly added parameters (not originally in JSBACH: τ, q, g₀, g₁) were derived from a synthesis of literature values. Most of the parameter ranges (limiting values for the parameters) were adapted from our previous work on a similar topic (Mäkelä et al., 2016). The parameter grouping was done to enhance optimisation, and the mechanism is explained in Sect. 2.7. Group I consists of parameters most directly affecting photosynthesis, group II parameters are intimately involved with soil moisture, and group III are the LoGro-P model parameters. The equations governed by these parameters are presented in Appendix A.

(Louis, 1979)

Table 2 Descriptions of model parameters with default values, range of acceptable values and references to equations in the paper or in the appendices. Parameters in the same group were calibrated simultaneously.

The start of the growing season in the JSBACH model is defined by a “spring event” in LoGro-P (Appendix A3) that induces leaf growth. The phenology model calculates a sum of ambient temperature (heatsum) since last autumn that is above the cutoff value T_alt, presented in Eq. (A10). It also calculates a variable threshold, defined in Eq. (A12), for the heatsum to reach. The threshold decreases based on the number of days the ambient temperature is below T_alt, whereas the heatsum increases. When the heatsum reaches the threshold, the plant leaves are free to grow.

However, coniferous evergreen trees do not shed all of their leaves for winter and the existing foliage enables them to quickly initiate photosynthesis in the following spring. The start of the photosynthetically active season in the model has been observed to occur too early in the boreal region by, e.g., Böttcher et al. (2016). In order to correct this behaviour, i.e. to restrain the respiration and photosynthesis of conifers in the early spring, we utilise a delayed effect of temperature for photosynthetic activity, introduced by Mäkelä et al. (2004). To calculate the reduction, we must first define the state of photosynthetic acclimation that Mäkelä et al. (2004, p. 371) present as “an aggregated measure of the state of those physiological processes of the leaves that determine the current photosynthetic capacity at any moment”.

The state of acclimation (S) is calculated from air temperature (T) with a delay prescribed by parameter τ (this is similar to the calculation of T_S in Appendix A14). S is then inserted into sigmoidal relation Eq. (1) to calculate a factor γ, a formulation that is adapted here from Kolari et al. (2007). Finally, γ is used to reduce the photosynthetic efficiency in Eq. (A1). T_1∕2 denotes the inflection point where γ reaches half of γ_max, k is the curvature of the function and γ=1 when S≥10.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}={\displaystyle \frac{T-S}{\mathit{\tau }}},\mathit{\gamma }={\displaystyle \frac{{\mathit{\gamma }}_{\max }}{\mathrm{1}+e^{k(S-T_{\mathrm{1}/\mathrm{2}})}}}\end{array}$$

The JSBACH model was also modified to include altogether six different stomatal conductance formulations following Knauer et al. (2015). These formulations include the pre-existing Baseline and Bethy versions as well as the Ball–Berry model and three of its variants. Model information is gathered in Table 3 for easy referencing and the detailed formulations are given in Appendix B. The limits of the slope of the stomatal conductance formulation parameter (g₁) were set to reflect commonly observed values from physiological measurements (Egea et al., 2011). The limits of $g_{\mathrm{1}}^{\mathrm{USO}}$ reflect the results presented by Lin et al. (2015).

We have also included two additional parameters (a and d in Table 2) for the Friend and Kiang (Friend and Kiang, 2005) stomatal conductance formulation in Eq. (B3). These parameters were not originally included in the optimisation, but the resulting cost function (Eq. 9) values were poor when compared to the other formulations. At that point, these parameters were included in the optimisation process. This increases the degrees of freedom for the Friend and Kiang model by two and therefore may give it an advantage when compared to the other Ball–Berry type formulations, which has to be considered in the interpretation of the results.

All of the stomatal conductance models contain an empirical water stress factor β, which reduces stomatal conductance as a function of volumetric soil water content (θ).

$$\begin{array}{}\text{(2)}& \mathit{\beta }=\left\{\begin{array}{ll}\mathrm{1},& \mathit{\theta }\ge {\mathit{\theta }}_{\mathrm{tsp}}\\ \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{pwp}}}{{\mathit{\theta }}_{\mathrm{tsp}}-{\mathit{\theta }}_{\mathrm{pwp}}},& {\mathit{\theta }}_{\mathrm{pwp}}<\mathit{\theta }<{\mathit{\theta }}_{\mathrm{tsp}}\\ \mathrm{0},& \mathit{\theta }\le {\mathit{\theta }}_{\mathrm{pwp}}\end{array}\right.\end{array}$$

In JSBACH, the stomatal conductance (g_s) is primarily resolved to estimate carbon fixation. The same g_s is then later used to calculate transpiration (Eq. A8). In the original JSBACH formulation (i.e. the Baseline version), the g_s is first resolved for unstressed canopy and then scaled by the water stress factor β. The Bethy approach is similar, but the conductance can also be limited by water supply (Eq. B2). In cases when the water supply is not the limiting factor, the calculations are similar to the Baseline version. In all of the empirical Ball–Berry variants, the stomatal conductance can be written as $g_{\mathrm{s}}=g_{\mathrm{0}}+c\mathit{\beta }{g}_{\mathrm{1}}$. The residual conductance (g₀) and the slope of the function (g₁) are both formulation-specific parameters as well as the factor c, which incorporates net photosynthesis and effects of atmospheric humidity and CO₂ concentration. The parameters g₀ and g₁ are part of our sampling and optimisation processes (group I in Table 2 when applicable).

The water stress factor (β) limits the carbon fixation and transpiration via the stomatal conductance formulation. Following Egea et al. (2011), it is also used to directly limit the net assimilation rate (A_n), as seen in Eq. (A1). The additional scaling (or limiting) factor for A_n takes the form β^q, so it is a function of both soil water content θ and the parameter q. Maximal reduction is achieved when q=1, and the reduction factor reverts to β. The minimal reduction occurs when q=0 and the reduction factor resembles a step function (at θ=θ_pwp). For any other value of q, it is a continuous convex function between the two extremes ${\mathit{\beta }}^q:[{\mathit{\theta }}_{\mathrm{pwp}},{\mathit{\theta }}_{\mathrm{tsp}}]\to [\mathrm{0},\mathrm{1}]$.

Knorr (1997)Knorr (2000)Ball et al. (1987)Leuning (1995)Friend and Kiang (2005)Medlyn et al. (2011)

Table 3 Stomatal conductance models with default values and range for g₁ and references to equations in Appendix B as well as related articles.

The ^∗ symbol indicates the Ball–Berry model and its variants.

2.4 Model simulations

The site-level measurements, used as model inputs, are air temperature, air pressure, precipitation, humidity, wind speed and CO₂ concentration as well as short- and longwave and potential shortwave radiation. Additionally, ET and GPP, derived from the EC measurements, are used to constrain and evaluate the model (as explained later in Sect. 2.8 and 2.9). We drive the model with half-hourly data but output daily values.

The initial state of the JSBACH model can be generated from predefined values of state variables (usually empty initial storage pools) or the model can be restarted from a file describing the state of some previous run. Depending on the area of interest, a model spin-up may be required to bring the model into a steady state. In our simulations, some of the more slowly changing variables (e.g. soil water content and LAI) need to be equilibrated, so a spin-up is required. This can be achieved by running the model over a set of measurements multiple times, each time restarting from the final state of the previous run.

The calibration period consists of the first 5 years given for the calibration sites in Table 1. The spin-up is achieved by looping over these 5 years, altogether four times (20-year spin-up), and then saving the state of the model at the end of the run. The actual calibration is started from the beginning of the calibration period, using the previously saved state variables. To reduce any bias this induces, the first year in the calibration run is removed from the cost function calculations. The spin-ups for the validation sites in Table 1 are similarly generated.

During the summer 2006, the Hyytiälä (FI-Hyy) measurement site suffered from a severe drought (Gao et al., 2017), leading to visible discolouration of needles. These events are difficult for models to capture and hence are of interest to modellers. We have previously and unsuccessfully attempted to optimise the JSBACH model (Mäkelä et al., 2016) for this event. Here we focus directly on the extended dry period (190–260th day of the year in 2006), during which the actual drought is mostly in effect between 210 and 235th DOY (day of the year). We adjusted some of the parameter values as those uncovered by the more general calibration, presented above. The spin-up was the same as for the calibration period, but at the end of the spin-up, the model was run forward to the start of the year 2006. Only values between the 190 and 260th DOY in 2006 were used in constraining the model.

2.5 Sampling process

We describe the modelling setup with the equation $\mathbf{y}=\mathcal{M}(\mathit{\theta },\mathbf{x})+\mathbf{e}$, where the aim is to reproduce the observations (y) with our model (ℳ), the driving data (x) and the current parameter values (θ). The residuals (e) depict how well the model reproduces the observations, and they form the basis of the likelihood function (formulated in Sect. 2.9), which is used to derive the parameter posterior distributions.

Using Bayes' rule on conditional probability we can write the parameter posterior density (p(θ,ℳ|x)) as a function of the likelihood (ℒ(x|θ,ℳ)), parameter prior distributions (π(θ)) and the model evidence (Z(x|ℳ)). As usual and from here on, we do not write ℳ in the Bayes' formula:

$$\begin{array}{}\text{(3)}& p\left(\mathit{\theta }\mathrm{|}\mathbf{x}\right)={\displaystyle \frac{\mathcal{L}\left(\mathbf{x}\mathrm{|}\mathit{\theta }\right)\mathit{\pi }\left(\mathit{\theta }\right)}{Z\left(\mathbf{x}\right)}}.\end{array}$$

We can now utilise the posterior density as a probability density function (pdf) for the parameters and infer the expectation values:

$$\begin{array}{}\text{(4)}& E\left[{\mathit{\theta }}_{\mathit{i}}\right]={\displaystyle \frac{\mathrm{1}}{Z}}\int {\mathit{\theta }}_{\mathit{i}}p\left(\mathit{\theta }\mathrm{|}\mathbf{x}\right)\mathrm{d}\mathit{\theta },Z=\int p\left(\mathit{\theta }\mathrm{|}\mathbf{x}\right)\mathrm{d}\mathit{\theta }.\end{array}$$

Above θ_i is the ith element of the parameter vector. Generally, Eq. (4) cannot be analytically solved; hence it is usually estimated numerically. Commonly this is achieved by one of the many MCMC methods, but in this study we apply APIS defined by Martino et al. (2015). APIS (Martino et al., 2015) is a Monte Carlo (MC) method that utilises a population of importance samplers (IS) to jointly estimate the target pdf (p(θ|x)) and the normalising constant (Z(x)) by a deterministic mixture approach (Veach and Guibas, 1995; Owen and Yi, 2000), whereas the MCMC methods do not care about the value of Z. We denote the importance sampling density as q(θ).

$$\begin{array}{}\text{(5)}& E\left[{\mathit{\theta }}_{\mathit{i}}\right]={\displaystyle \frac{\mathrm{1}}{Z}}\int {\mathit{\theta }}_{\mathit{i}}r\left(\mathit{\theta }\right)q\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta },\text{where}r\left(\mathit{\theta }\right)={\displaystyle \frac{p\left(\mathit{\theta }\mathrm{|}\mathbf{x}\right)}{q\left(\mathit{\theta }\right)}}\end{array}$$

Above r is the reweighing factor that is the driving force in importance sampling. We will next give a summary description of the sampling process with comparison to a general multichain MCMC approach (since MCMC methods are more commonly used in these types of situations).

1. The initialisation of a multichain MCMC sampler and APIS are very similar. In our simulations, APIS is set up as 40 simultaneous and independent importance samplers. This is similar to an independent 40-chain MCMC sampler. Each sampler or chain has a random starting location drawn from a uniform distribution defined by the parameter ranges, given in Table 2. The initial sampling (or prior) distribution for each sampler is also randomly generated – we use truncated Gaussian distributions with diagonal covariance matrices, where the standard deviations are randomised. The sampling distributions will evolve throughout the process.

2. In an MCMC setup, the model would be run once (for each chain) and evaluated and then the draw (parameter values) accepted or rejected accordingly. In APIS, instead of a single element (one run) we use a sample size of 50. This means that we draw 50 elements with each IS sampler (or “chain”) independently. These draws are then evaluated and reweighted as presented in Eq. (5).

3. The 50 reweighted draws (for each IS sampler separately) are used to calculate a new location for the sampling distribution. This location is automatically accepted (no rejection criteria), and we also adapt the shape of the distribution using the self-normalising estimator by Cornuet et al. (2012).

4. Additionally, all of the draws in APIS are used to calculate “global” estimates of the parameter expected values. This process utilises the deterministic mixture approach (Veach and Guibas, 1995; Owen and Yi, 2000) and is fully iterative with no need for any recalculations as the previous estimates are directly adjusted (no information is lost either).

MCMC chains track the evolution of single elements and occasionally adjust the sampling distribution. The sample size in APIS is larger (it is not a Markov chain method) and the focus is on the evolution of the locations of the sampling distributions, not on the individually drawn elements. These location parameters are expected to be around all the modes of the target and the deterministic mixture ensures the stability of the estimation of the (global) parameter expected values. As an importance sampler, APIS is also a variance reducing method.

Before taking a more detailed look at APIS, we make some further notes about the sampling process. The first element of the 50 draws (item 2 in the list above) is always fixed as the current mean. We run the spin-up (Sect. 2.4) and generate the model starting state only for the proposal means and use the same state for the other 49 draws (perturbed around the proposal mean). This requirement stems from a need to reduce computational time as running the model to a steady state is costly. This approach might induce some discrepancies, but they are mitigated by removing the first year of the calibration simulations (as explained in Sect. 2.4). We also slightly reduce the importance weights of the 49 samples (more reduction for samples further from the proposal mean), when calculating the new location parameters (item 3 in the list above) – the reduction only (slightly) slows the adaptation of the IS sampler locations. Finally, we note that this approach ensures that we run the proposal means, which are the focus in APIS, with the correct spin-up.

2.6 Adaptive population importance sampler

Normally, only the location parameters of the IS proposals are adapted, but we also adapt the shape parameters using the self-normalising estimators by Cornuet et al. (2012). APIS is able to utilise different or a mixture of normalised proposals densities, but we use truncated Gaussian proposals with diagonal covariance matrices.

In our simulations, APIS is formed of 40 independent IS estimators. Each estimator draws a sample ${\mathit{\theta }}_i,i\in \mathit{\{}\mathrm{1},\mathrm{\dots },N\mathit{\}}$, of size N=50 at a time from their own proposal distribution $q_j\left(\mathit{\theta }\right),j\in \mathit{\{}\mathrm{1},\mathrm{\dots },M\mathit{\}},M=\mathrm{40}$. The estimator then calculates the importance weights ($w_{ij}=\frac{p\left({\mathit{\theta }}_i\mathrm{|}\mathbf{x}\right)}{q_j\left({\mathit{\theta }}_i\right)}$) for each sample. The location (μ_j) and shape (C_j) parameters (Cornuet et al., 2012) of each proposal are updated using only samples (and weights) drawn from q_j. The new shape parameters are formed as a mean of the previous estimate and C_j, as calculated below.

$$\begin{array}{}\text{(6)}& {\mathit{\mu }}_j={\displaystyle \frac{{\sum }_i{w}_{ij}{\mathit{\theta }}_i}{{\sum }_i{w}_{ij}}},{\mathbf{C}}_j={\displaystyle \frac{{\sum }_i{w}_{ij}({\mathit{\theta }}_i-{\mathit{\mu }}_j)({\mathit{\theta }}_i-{\mathit{\mu }}_j{)}^T}{{\sum }_i{w}_{ij}}}\end{array}$$

The simple IS estimators alone are rarely sufficient if the target is even slightly complicated. One classical way of tackling this problem is to join multiple IS estimators together. The simplest approach is to calculate the weights for each of these estimators separately and to normalise the result by the combined sum of all weights. However, this leaves the estimators susceptible to “bad” proposals. APIS suppresses the bad proposals by utilising the deterministic mixture approach (Veach and Guibas, 1995; Owen and Yi, 2000) presented in Eq. (7), where each proposal q_j is evaluated at all the drawn samples and weighted by the amount of samples drawn (N_j=50) from that proposal. This is equivalent to joining the normalised proposal densities together and evaluating the joint pdf.

$$\begin{array}{}\text{(7)}& w_{ij}={\displaystyle \frac{p\left({\mathit{\theta }}_{ij}\mathrm{|}\mathbf{x}\right)}{{\sum }_j\left(\frac{N_j}{{\sum }_k{N}_k}\right)q_j\left({\mathit{\theta }}_{ij}\right)}}\end{array}$$

The parameter expectation values and the normalising constant in Eq. (5) can now be estimated by Monte Carlo integration using weights calculated in Eq. (7).

2.7 Parameter optimisation

The APIS algorithm is a rather robust method meant for examining the full target probability distribution and locating the modes of the target distribution. Adaptation in APIS utilises multiple draws simultaneously, which can easily lead to few parameters controlling this process (the marginal density of one or few parameters dominates the calculations). Since we also did not run the model spin-up for all drawn samples (although the discrepancies should be minimal), we utilise a simple custom stochastic optimiser to locate the optimal set of parameter values. This optimiser is run after the APIS calibration simulations and separately for the drought period. The optimiser utilises the exact same datasets (calibration, validation, observations, etc.) as APIS, the spin-up is generated for all drawn samples separately and the initial state of the algorithm is the mean value of the APIS final configuration (location parameters).

Our optimiser is a simple random sampler amplified by the “velocity” of the last jump (the idea is similar to Hamiltonian or Hybrid Monte Carlo by Duane et al., 1987). We draw a set of samples from a small Gaussian proposal distribution in the vicinity of the current best estimate and calculate the cost function for the samples. Whenever a better point is found (smaller cost function), we jump to that (update the mean of the proposal distribution). The velocity of the jump (for us merely the distance of change in each parameter) is then added to the new mean (with a maximal limit of 1 standard deviation in the proposal distribution), but it is reduced and eventually removed if a better sample is not found.

The covariance matrix of the proposal distribution is recalculated at predefined intervals (for all parameters). Additionally, we utilise a subset sampling procedure, where the samples are first drawn from the full parameter space; in the next step they are drawn only from group I in Table 2 (the rest are kept at their current optimal values), followed by groups II and III and then back to the full parameter space. When the number of parameters is reduced, we are more likely to find a better set of parameter values. We have kept the parameters mostly affecting the same processes in the same group, but some dependencies may not be apparent, and hence it is also important to draw samples from the full parameter space.

Figure 1 Examples of the evolution of the APIS algorithm from the Bethy calibration. The left panel is the kernel density estimate of the location parameters at the start of the process (black), after 20 iterations (blue) and after 100 iterations (green). The right panel shows the location parameters (grey), their mean (red) and 1 standard deviation (dashed) as well as the global estimate (yellow, calculated with the deterministic mixture approach) of the parameter expected value.

2.8 Simulation analysis

Even though APIS is not a Markov chain method, we can (naively) interpret the evolution of the location parameters of each IS sampler as chains. The resulting 40 chains have random starting positions, but they are relatively short (we present results from the Bethy calibration, where the chains were adjusted 100 times); hence we did not discard any of the samples. We test the convergence of these chains with the Gelman–Rubin diagnostic tests (Gelman and Rubin, 1992), comparing the variance between the chains to the variance within each chain and calculating the potential scale reduction factors ($\widehat{R}$). We also test the stability of the (parameter) global expected value estimate (using the deterministic mixture approach) by calculating the difference of the final global expected value and the mean of the location parameters (at each iteration). We denote this test as δ and report the number of the iterations when this difference is below 5 % of the parameters range, given in Table 2.

In order to visualise the results, we have utilised a Gaussian kernel density estimation (KDE) to produce distributions from the APIS simulation location parameters. In practice, KDE places a Gaussian distribution centred at each sample, and the constructed composite distribution is an estimate of the underlying actual distribution. The bandwidth for the distributions is calculated using Scott's rule (Scott, 2004): the data covariance matrix is multiplied by a factor $n^{\frac{-\mathrm{1}}{d+\mathrm{4}}}$, where n is the number of data points and d is the number of dimensions.

The effectiveness of each parameter was calculated from the final state of each optimisation process. This was done by first setting all parameters to their optimised values. Then we (evenly) sampled each parameter separately from their range of acceptable values, given in Table 2, and calculated the corresponding cost functions. For each parameter the maximum difference in these cost function values (and the optimised value) was recorded. The parameters (within each optimisation) were then ordered by these numbers (with the highest difference meaning highest effectiveness) and separated into three groups: highest (most effective) and lowest (least effective) effectiveness values and the rest. This effectiveness relates to how the APIS “sees” the sampling process – the 50 draws are evaluated simultaneously, and a very effective parameter can easily mask the influence of a less effective one (the marginal density of one or few parameters dominates the calculations).

We report the slope of the regression line (b) and the coefficient of determination (r²) between the observations (y_i) and the model output (x_i). The slope of the regression line is highly indicative of the model bias (the difference of the expected values of the observations and the model). Hence we interpret the bias directly from b (in our results the regression lines pass near origin so the differences this induces are negligible).

$$\begin{array}{}\text{(8)}& b={\displaystyle \frac{{\sum }_i(x_i-\stackrel{\mathrm{‾}}{x_i})(y_i-\stackrel{\mathrm{‾}}{y_i})}{{\sum }_i(y_i-\stackrel{\mathrm{‾}}{y_i}{)}^{\mathrm{2}}}},r^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{{\sum }_i(x_i-y_i{)}^{\mathrm{2}}}{{\sum }_i(y_i-\stackrel{\mathrm{‾}}{y_i}{)}^{\mathrm{2}}}}\end{array}$$

2.9 Cost function

The Bayesian framework requires a likelihood function that optimally combines pointwise model and observational errors. The JSBACH model error is unknown as is the (pointwise) observation error. We could use a general type of error estimate (such as that of 20 % of the flux value) for the observations but would have to include a minimal site and instrumentation-dependent precision. In this study, the full error is treated as Gaussian white noise. Because of these limitations, we are calling and defining our likelihood as a cost function. It is calculated with the same parameter values for each site, using site-specific daily measurements with the gap-filled, low-quality and winter (between the 315th and the 75th day of the year) values removed (resulting in N_ET and N_GPP points). These site-level estimates are averaged to produce the actual cost function, which is then returned for the algorithm to produce an estimate that is independent of the characteristics of any single site.

The cost function (Eq. 9) in our simulations is based on the normalised mean squared error (NMSE) estimates of the daily GPP and the daily ET. The residual of each variable is divided by the mean of observations, as has been previously done by, e.g., Mäkelä et al. (2016), Knauer et al. (2015), Groenendijk et al. (2010) and Trudinger et al. (2007). We make use of this approach since we needed to balance two series of different magnitudes (ET and GPP). The residuals are additionally divided by the (site-specific) number of observations so that the cost function is not biased towards any-specific site. The cost function (without the normalisation) can be interpreted as a negative log-likelihood function with a (Gaussian) error term equal to the observational mean.

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mathrm{cf}}_{\mathrm{1}}& =\stackrel{{\text{NMSE}}_{\mathrm{ET}}}{\overbrace{{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{ET}}}}\sum {\left({\displaystyle \frac{{\mathrm{ET}}_{\mathrm{mod}}-{\mathrm{ET}}_{\mathrm{obs}}}{\stackrel{\mathrm{‾}}{{\mathrm{ET}}_{\mathrm{obs}}}}}\right)}^{\mathrm{2}}}}\\ & +\stackrel{{\text{NMSE}}_{\mathrm{GPP}}}{\overbrace{{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{GPP}}}}\sum {\left({\displaystyle \frac{{\mathrm{GPP}}_{\mathrm{mod}}-{\mathrm{GPP}}_{\mathrm{obs}}}{\stackrel{\mathrm{‾}}{{\mathrm{GPP}}_{\mathrm{obs}}}}}\right)}^{\mathrm{2}}}}\end{array}\end{array}$$

We also use a modified version of this cost function, where the NMSEs are weighted by factors based on coefficients of determination (r²) defined in Eq. (8). This latter cost function is only used during the separate drought period optimisation for Hyytiälä. During the drought we are more interested in the correct timing of the change in GPP and ET fluxes, rather than the size of the actual change. The aim is to correctly reproduce the changes in the water use efficiency (WUE) of plants, which we interpret here as the pointwise ratio of (ecosystem level) GPP to ET. The NMSE values ensure that the overall amplitude of the fluxes will remain satisfactory.

$$\begin{array}{}\text{(10)}& {\mathrm{cf}}_{\mathrm{2}}=(\mathrm{1}-r_{\mathrm{ET}}^{\mathrm{2}}){\text{NMSE}}_{\mathrm{ET}}+(\mathrm{1}-r_{\mathrm{GPP}}^{\mathrm{2}}){\text{NMSE}}_{\mathrm{GPP}}\end{array}$$

3 Results

First we present the performance of the APIS algorithm and the parameters themselves, followed by site and stomatal conductance model-specific results and finally an examination of the Hyytiälä drought event in 2006. For simplicity, we use the name of the stomatal conductance model to refer to the JSBACH model utilising that stomatal conductance formulation.

The evolution of the APIS algorithmic process is presented in Fig. 1 for three parameters from the calibration of the Bethy model. The chosen parameters highlight different levels of identifiability for the algorithm (with the given cost function). The first parameter (f_C3) shows a well-identifiable situation, where the algorithm quickly locates the area of high probability. The second parameter (θ_dr) is also identifiable, but the speed of convergence is diminished. The last example (C_decay) represents situations where the parameter is not constrained. We have included images of the APIS chains for the other parameters as Supplement S1 along with parameter posterior estimates at 20 iterations with the Bethy and Ball–Berry formulations.

We also report the results of the Gelman–Rubin (Gelman and Rubin, 1992) and δ tests in Table 4. Both of these tests indicate that the algorithm is performing well at 20 iterations – the values of $\widehat{R}\approx \mathrm{1}$, which means that further simulations are unlikely to improve the variance estimates. However, for some parameters, the convergence of the global estimate is slow (as also seen in the Supplement Fig. S1 for, e.g., τ, c_b and q). The APIS sampling process did not reveal any multimodal distributions and thus provided suitable initial conditions for the optimisation.

Table 4 Parameter scale reduction $\widehat{R}$ (at APIS iteration) and stability δ (threshold number of iterations) estimates from the Bethy simulations.

3.1 Optimised parameters

The results of the optimisation process are gathered in Table 5. There is an overall agreement on the values of the most prevalent parameters (see the bold and italic characters in Table 5 between the models). Most notably, the permanent wilting point (θ_pwp) and the point above which transpiration is unaffected by soil moisture stress (θ_tsp) have been significantly lowered. The LoGro-P parameters, which affect the timing of the spring and autumn events, are expected to contribute only little to the cost function. The coniferous evergreen trees do not shed all their leaves for winter, and therefore the timing of the bud burst is not as critical as for, e.g., deciduous trees. Additionally, because of the existing foliage, the state of acclimation parameter τ that depicts the reduction in carbon assimilation in the early spring likely dominates the phenology parameters that determine when new leaves start to grow.

Table 5 Parameter default and optimised values for the calibration period with corresponding cost function value. The values written in boldface were the most effective and the italic values the least effective for the given experiment. Also presented are the fixed parameter values for the drought period optimisation, with “opt” referring to the use of the corresponding optimised value from this table.

Some of the parameters have converged to their limiting values, which can reflect deficiencies in the model structure or the preset parameter ranges. Convergence to the boundary can also be a problem in model calibration, but in this experiment, the algorithms were able to cope with the situation as APIS located the area of high probability and the optimiser located the maxima. The different parameter effectiveness levels reported in Table 5 can be roughly equated to the identifiability situations in Fig. 1. The effectiveness levels are highly situational (e.g. they depend on the sampling limits in Table 2 given for each parameter) and merely reflect the parameter identifiability in the APIS process. Low effectiveness complements the test results in Table 4, as the tests may indicate good performance for a parameter (e.g. for S_range) that is ineffective in the simulations.

3.2 Annual cycles

We present the average annual cycles for the validation period and for all sites in Fig. 2 using the Bethy formulation that is part of the standard model. The annual cycles generated with the other stomatal conductance models are added as Supplement S2. The parameters of the regression lines (b and r²) between the measured and modelled ET and GPP fluxes of all the models are gathered in Table 6. These indicators have been calculated using all corresponding values regardless of the quality of the data. The sites are in the same order as in Table 1 with the six calibration sites first, followed by the four sites used only for validation. We have also included a supporting synthesis of the b and r² values between the model simulations with the default and optimised parameter values as Supplement S3.

Figure 2 Validation period average annual cycles of evapotranspiration and gross primary production; observations (black) and the model using the Bethy stomatal conductance formulation with default (green) and optimised (blue) parameterisation. Also presented are daily model values cross plotted against observations with corresponding slope of the regression line (b) and the coefficient of determination (r²).

Table 6 Slope of the regression line (b) and the coefficient of determination (r²) for the different stomatal conductance formulations during the validation period with the optimised parameters. We have written the best values of b and r² in boldface for each site and italicised the abbreviations of the separate validation sites.

The optimisation has improved the model bias and the correlation coefficients for the GPP in Fig. 2 for nearly every site, with the exception of deteriorating bias for Poker Flat (US-Prr) and Zotino (RU-Zot). Additionally, the improvement in the timing of the springtime increase in the GPP is apparent. All of the correlation coefficients for the ET in Fig. 2 have also been improved but the model bias has mostly increased.

3.3 Drought event

The resulting parameter values, from the optimisation during the drought conditions in Hyytiälä (FI-Hyy) in the summer of 2006, are presented in Table 7. Setting the maximum carboxylation rate to a constant value ($V_{\mathrm{C},\max }=\mathrm{52.0}$) enabled the full use of the dynamical range of q – the idea was to ensure that V_C,max does not dominate the optimisation, any value for q is possible, and it is able to influence the outcome. The LoGro-P parameters and τ were fixed to their optimised values, presented in Table 5, as they should not be affected by the drought. Likewise, the values of other parameters (not presented in Table 7) were set as compromises between the stomatal conductance formulations.

We can now compare the parameter values in Table 7 to those in Table 5. The values of the relative humidity parameter (θ_hum) and the residual stomatal conductance (g₀) have remained nearly unchanged, but the rest of the parameters have quite varied values. The leaf internal-to-external CO₂ concentration (f_C3) as well as the slope of the stomatal conductance (g₁) are at the lower bound (expect g₁ for BB – Ball–Berry). Noticeably, the unified stomatal optimisation model (USO) only changes the values of θ_tsp and q and leaves the rest of the parameters almost untouched.

The changes these different parameterisations have on the model output are visualised in Fig. 3. All of the stomatal conductance models, with default parameterisation, suffer from ET values that are too low before (and during) the actual drought. This behaviour was corrected during the general optimisation but has partially re-emerged with the dry-period parameters for the Baseline, Ball–Berry, Leuning, and to a lesser degree the Friend and Kiang formulations. Most of the models also exhibit ET values that are too high during the actual drought with the generally optimised parameter values. This behaviour was also corrected with the dry-period optimisation, but the Baseline and especially the Bethy model now suffer from too strong a drawdown of ET. These models also demonstrate the drawdown that is too strong for the GPP. The GPP itself was greatly improved with both optimisations and for all models. The dry-period optimisation of the USO also managed to correct the erroneous GPP of the general optimisation during the actual drought, whereas the GPP of other formulations has remained roughly the same as with the general optimisation. The USO formulation results in the best fits for r² and b with the dry-period optimisation.

The Bethy and the USO models demonstrate the most variability in the β-function values in Fig. 3 (rightmost panels), for the dry-period optimisation. We selected these two stomatal conductance formulations to examine the changes to the WUE of plants during the extended dry period. The highlighted observations in Fig. 4c and f show a clear path of development for the drought where the observations imitate the letter δ. The colourings follow the β-function values in Fig. 3 between the red vertical lines. Both observational colourings (same as the model colouring) are similar and depict, initially, a linear decrease in both ET and GPP, followed by a rapid decline in ET and a delayed decline in GPP. The recovery of plants from the drought can also be seen as the colouring starts to turn lighter. The models depict a more linear response of GPP to ET as the drought develops, although with USO we can see a few more similarities in the pattern of the values.

Table 7 Optimised parameter and corresponding cost function values with different stomatal conductance formulations for the extended dry period.

Figure 3 Hyytiälä site drought in summer 2006. The time series for evapotranspiration and gross primary production are 5 d running averages and for β-function daily values. The observations are plotted in black and the model with default parameterisation in green, calibration period optimisation in blue and the dry-year optimisation in magenta. The red vertical lines indicate the start and end of the actual drought.

Figure 4 Hyytiälä site water use efficiency for the Bethy and USO formulations. Scatter-plotted are the dry-period 5 d running averages of ET and GPP, coloured by the intensity of the drought (β function). Panels (a, d) depict the model with the more generally optimised parameter values and (b, e) with the drought optimisation, and (c, f) present the corresponding observations, coloured by the same intensity as in (b, e). The grey points are from the corresponding time during the 2 previous years.

Finally, we used both optimised parameter sets (Tables 5 and 7) to produce the ET and GPP cycles for all sites and stomatal conductance models. This analysis (not shown) verified that in general conditions, the Table 5 parameter values produced better estimates in general. The b and r² values for the ET were systematically better for all stomatal conductance formulations (except one). There was some variation in the indicators for the GPP, where approximately a third of the values (of mostly r²) are better with the dry-period parameter set. These differences are mostly attributed to increased model bias (decreased b) that is explained by the lower values of g₁. Overall, the more general optimisation provided systematically better or comparable results to the dry-period optimisation. The exception is the USO formulation, which had an approximately 1:1 distribution of best values for both variables in between the parameter sets.

4 Discussion

We will first discuss the validity of our approach and the simulation setup, followed by an examination of the success of the modifications made to the model, and close with some further remarks on the parameter values.

4.1 Validity of the simulations

Before we calibrated the model, we fixed the limiting value for LAI and adjusted the site-specific vegetative area fractions to reproduce the measured site-level maximum of LAI. In the simulations, we focused on boreal coniferous forests, where light penetration is deep and the light conditions are homogenous – consequently we could assume a homogenous leaf distribution. Furthermore, the JSBACH model takes into account leaf clumping, and we can assume the leaf orientation and shape to be similar throughout the study sites. Therefore, we argue that reproducing the site-level maximum of LAI is appropriate approach in this study. Together with parameter calibration, it has resulted in improved ET and GPP fluxes as can be verified from the b and r² values in Fig. 2. The improvements in b and r² are mostly seen in the GPP flux, which can be explained by the fact that the stomatal conductance in JSBACH is primarily resolved for carbon assimilation, and the same conductance is then used for transpiration (Eq. A8). Additionally, GPP is derived from the EC measurements by flux partitioning – this tends to remove some of the flux instabilities (that are still present in the ET).

We encountered difficulties in reproducing the fluxes for the validation sites with low LAI (i.e. RU-Zot and US-Prr). This can be a consequence of the area scaling as the adjustment linearly changes the proportions between vegetative area and bare soil. Another reason is the lack of the site understorey in these simulations. For example, approximately half of the CO₂ fluxes (and consequently roughly half of the GPP) for Poker Flat are produced by the site understorey (Ikawa et al., 2015). Additionally, there are also many parameters describing site-specific soil properties (such as porosity) that were not part of the optimisation and may be inaccurate. These effects may also be pronounced due to the changes in parameters affecting soil moisture as well as the area scaling.

There were no clear differences between sites dominated by pine or spruce. Neither did we notice any particular effect on the bias, NMSE or correlation coefficient that could be explained by geographical location, stand age,s or annual precipitation or temperature. We optimised the model for individual (calibration) sites as well (not shown). Mostly this changed the values of parameters (such as V_C,max and g₁), affecting the amplitude of the modelled fluxes. These parameters can be viewed to be more site-specific, a characteristics that is reduced in a multi-site calibration – the possibility of highly site-specific properties (and parameter values) can also explain the difficulties in reproducing the validation site observations. We are omitting these results as single-site optimisation can be viewed as overfitting the model and the results do not provide any additional insights.

The APIS performance tests (Gelman–Rubin and δ) indicate that the algorithm is performing well at 20 iterations, but the convergence of the global estimate for some parameters is slow. This is mostly a direct result of the normalisation of the cost function that inflates the target distribution, which reduces the parameter sensitivity to observations and gives too much weight to the initial locations and draws. Without the normalisation, the algorithm would also converge faster. Additionally, APIS is meant to examine the full target distribution with only some sequentiality – 20 iterations (or less) should be sufficient for APIS to locate the modes of the target. In longer APIS simulations, the global estimate would likely benefit from discarding the first half of the samples, but this would require the estimate to be recalculated at each iteration (from the drawn samples) as it could not be calculated iteratively.

4.2 Delayed effect of temperature

We modified the JSBACH model by introducing the delayed effect of temperature for photosynthesis to restrain the respiration and photosynthesis of conifers in spring. The effect of this (delayed increase in GPP) is apparent in the annual GPP cycles of CA-Qfo, FI-Hyy, FI-Ken, FI-Sod and RU-Zot in Fig. 2. The delay is in place for the other sites as well, but the effect is less apparent in the figure. This delay is to a lesser extent also reflected in transpiration, and consequently in ET, as can be seen, e.g., at FI-Hyy and FI-Sod – for other sites this effect is not clear. The correction in the ET values can lead to an increase in model bias as is the case with Sodankylä (FI-Sod), where the autumn values that are too low in the default model were previously compensated for by springtime values that too high (in the sense of annual ET). Correcting the springtime behaviour leads to an increase in bias, but this should not be viewed as a fault in the optimisation as the model was previously mitigating an erroneous behaviour (too low an autumn ET) with another (too high a springtime ET).

Mäkelä et al. (2004) used a linear dependency of photosynthetic efficiency to the state of acclimation and reported 13.75 d to be the best fit for the adjustment period length (τ). Kolari et al. (2007) utilised a sigmoidal relation and reported the value of 8 d but noted that the range of values resulting in a good fit was large (5–10.4 d). Linkosalo et al. (2014) came to a similar conclusion when they encountered a near-flat distribution for τ in the range of 1–12 d. In our simulations τ exhibits larger optimal values (nearly 15 d), which is most likely due to the model adapting to the multi-site calibration (as sites have different characteristics, a longer acclimation period accounts better for these variations).

4.3 Stomatal conductance models

We examined the model behaviour with six stomatal conductance formulations, and the resulting b and r² values are presented in Table 6. The best performance (bolded values) in simulated ET is achieved by the BB model for bias and the Bethy formulation for r². These two models also share the best performance in the GPP bias, whereas the best r² values for the GPP are demonstrated by the F&K (Friend and Kiang) model, followed by the USO formulation. Calculating the number of best values demonstrated by each model, we obtain that the best performance is shared by the Bethy (12) and F&K (12) formulations, followed by the BB (11) model. However, we note that some of the “best values” are only marginally better than comparable values. Additionally, we used two more parameters (a and d) for the F&K formulation than for the other Ball–Berry formulations. Likewise, we could have, for example, included the factor D₀ (Eq. B3) in the optimisation, which would have likely improved the performance of the Leuning model. Similarly to the results by Knauer et al. (2015), based on this (general) calibration, there is no clear single candidate for the best stomatal conductance formulation.

The model behaviour was also examined during the Hyytiälä drought of 2006. Some of the parameter values were kept fixed during these simulations; most of the fixed parameters should not affect the drought period calibration, but there are exceptions, such as the maximum carboxylation rate V_C,max. It can be argued that, e.g., both the parameters V_C,max and g₁ should decrease (Egea et al., 2011) during the drought, but we decided to fix V_C,max to get a better response for q. The best fit to the observations was achieved by the USO formulation, as seen in Fig. 3, with remarkably similar parameter values to the general optimisation. The USO was also able to (somewhat) replicate the “δ” shape of the drought in Fig. 4.

The stomatal conductance function ($g_{\mathrm{s}}=g_{\mathrm{0}}+c\mathit{\beta }{g}_{\mathrm{1}}$) also incorporates the soil water parameters θ_tsp and θ_pwp in the form of the β function as portrayed in Eq. (2). The changes in the values of these parameters (mostly g₁,θ_tsp and θ_pwp) are intertwined. During the drought, the decrease in the optimised values of g₁ is expected as the plants close their stomata to minimise the loss of water by transpiration (Egea et al., 2011; Zhou et al., 2013). The same effect is also achieved by increasing the values of θ_tsp and θ_pwp as this decreases the values of the β function. The higher values of g₁ during the more general optimisation are better reflected by Franks et al. (2018), whereas the lower values during the drought are more in accordance with physiological observations by Egea et al. (2011). Likewise, Lin et al. (2015) found higher values for g₁ (both boreal area and gymnosperm trees) using the USO.

In general, the site-level estimates of (g₀ and) g₁ are sensitive not only to the stomatal conductance formulation but also, e.g., to the structure of the underlying model and the value of other parameters, such as maximum carboxylation rate (V_C,max). (Wang, 1996, Table 1, Control) reported g₁=3.78, using a Leuning model similar to ours, where $(\mathrm{1}+D_{\mathrm{S}}/D_{\mathrm{0}})$ is replaced by D_S. Thum et al. (2007) approximated $g_{\mathrm{1}}^{\mathrm{BB}}$ to be 5 for Sodankylä while estimating the variation in the values of V_C,max and maximum rate of electron transport J_max. We would suggest that the limiting values θ_pwp and θ_tsp should be optimised or fixed before introducing additional tuning factors such as mesophyll conductance or scaling the β in multiple ways in the stomatal conductance formulations (Egea et al., 2011). Our simulation setup for q corresponds to the configuration 5 (C5) by Egea et al. (2011), with variables q=q_B and fixed value q_S=1.

4.4 Parameter values

Some of the parameters in this study have been calibrated before by, e.g., Kattge et al. (2009) and Knorr et al. (2010). Our approach differs from these as we required the model to reproduce the site-level maximum of LAI. In contrast, e.g., Knorr et al. (2010) found the structural limit for (all-sided) LAI to be 4.2, which is considerably lower than the measured LAI for many of the sites in Table 1. Our approach directly scales the vegetative area, so it also scales GPP and also the amount of rain available for plants (as rain is directed to bare soil and vegetative area). This means that the parameter values should not be directly compared without taking the different paradigms into account. However, our optimised V_C,max values are in between 62.5 reported by Kattge et al. (2009) and 29.3 by Knorr et al. (2010) and are in line with the yearly cycle presented by Ueyama et al. (2016).

The exponential scaling factor q in Eq. (A1) of the β function (Eq. 2) was revealed to be ineffective in our optimisation as indicated in Table 5. In our simulations, this situation arises as the effective range of the β function has been lowered by reducing θ_pwp and θ_tsp. The actual soil moisture is rarely below the fraction θ_tsp, so q is constrained with a very limited number of datapoints and thus only has a minimal effect on the fluxes and the cost function. Therefore, the values presented for q in Table 5 can be unreliable and even unrealistic. This situation is remedied in the drought period optimisation when the soil moisture is low. The resulting values for q in Table 7 have a wide dispersion, although they are mostly at the lower end. This signifies that the additional GPP reduction is mostly gradual, with a steep decrease near the permanent wilting point θ_pwp.

The values of soil water parameters are closely grouped in the optimisations except for the values of θ_pwp during the drought. This can occur due to a larger impact of the different stomatal conductance formulations on the accumulating soil water content than assumed – this can also be seen from the differences in the β-function values in Fig. 3. Furthermore, the values of θ_tsp and θ_pwp have been considerably lowered from their default values in both optimisations. This change can be perceived in at least two different ways. Either the boreal forests are not generally limited by soil moisture stress (except in the case of extreme drought) or the water retention capabilities of the soil (in the model) have been systematically overestimated. The latter seems unlikely, in the light of results by, e.g., Gao et al. (2016).

5 Conclusions

APIS is a recent method, capable of estimating complicated multidimensional probability distributions using a population of different proposal densities. The algorithm was able to produce reasonably stable estimates for most parameters quickly. Prior to calibrating the model, we adjusted the site-specific vegetative area fractions to reproduce the measured site-level maximum of LAI. This practical approach resulted in improved ET and GPP fluxes, although we encountered difficulties in replicating these for sites with low LAI. The model parameters were optimised simultaneously for all sites without any additional site-level tuning. The parameters that were most effective in the optimisation processes were consistent for all stomatal conductance formulations.

The introduction of the S function, to delay the start of the vegetation active season, has corrected the springtime increase in GPP for conifers throughout the sites used in this study. The parameters θ_tsp and θ_pwp, which set the range for the soil moisture stress function β, were both systematically lowered and optimised to nearly identical values for all stomatal conductance models. The low effective range for the β function rendered the experimental parameter q nearly ineffective in the more general optimisation. The dry-period optimisation increased the effective range of the β function and the importance of q, which resulted in a highly nonlinear (additional) reduction in the net assimilation rate. Overall, this fact and both optimisations indicate that boreal forest transpiration is not limited by soil moisture stress under normal conditions.

The optimisation improved the predictive skill of the model with all stomatal conductance formulations as was seen during the validation period. The Bethy, Ball–Berry, and Friend and Kiang versions were the most in agreement with the observations, although the differences between these and the other formulations were small. Most of the model versions had some problems during the extended dry period, and the best b and r² values were achieved by the unified stomatal optimisation model. Additionally, the optimised parameter values of the USO for the dry period were the most similar (of all stomatal conductance formulations) to those of the more general optimisation.

Appendix A: Parametric equations within JSBACH

In this Appendix we present the most relevant equations that are governed by the parameters in Table 2. The Appendix is divided into sections that coincide with the parameter groups.

A1 Photosynthesis

The Farquhar model (Farquhar et al., 1980) is based on the observation that the assimilation rate in the chloroplast is limited either by the carboxylation rate (V_C), induced by the Rubisco enzyme, or the light-limited assimilation rate (J_E). The total rate of carbon fixation is reduced by the amount of dark respiration (R_d), resulting in net assimilation rate (A_n). The experimental scaling factor β^q (Egea et al., 2011) is based on soil moisture stress in Eq. (2), which takes effect (β<1) when soil moisture is significantly reduced. This scaling is used by all stomatal conductance formulations. We have also introduced here in equation form the actual reduction to photosynthesis by γ from the delay in the start of the vegetation active season in Eq. (1).

$$\begin{array}{}\text{(A1)}& A_{\mathrm{n}}={\mathit{\beta }}^q(min(\mathit{\gamma }{V}_{\mathrm{C}},J_{\mathrm{E}})-\mathit{\gamma }{R}_{\mathrm{d}})\end{array}$$

Oxygenation of the Rubisco molecule reduces the carboxylation rate, which is given as

$$\begin{array}{}\text{(A2)}& V_{\mathrm{C}}=V_{\mathrm{C},\max }{\displaystyle \frac{C_{\mathrm{i}}-{\mathrm{\Gamma }}_{\star }}{C_{\mathrm{i}}+K_{\mathrm{C}}(\mathrm{1}+O_{\mathrm{i}}/K_{\mathrm{O}})}}.\end{array}$$

Here C_i and O_i are the leaf-internal CO₂ and O₂ concentrations, Γ_⋆ is the photorespiratory CO₂ compensation point, and K_C and K_O are Michaelis–Menten constants parameterising the dependence on CO₂ and O₂ concentrations. Furthermore, leaf-internal CO₂ concentration depends on the external (ambient) concentration C_a (in the Baseline and Bethy formulations and unstressed conditions) by

$$\begin{array}{}\text{(A3)}& C_{\mathrm{i}}=f_{\mathrm{C}\mathrm{3}}{C}_{\mathrm{a}}.\end{array}$$

Likewise, the light-limited assimilation rate can be expressed as a function on electron transport rate (J), which is a function of radiation intensity (I) in the photosynthetically active band, the maximum electron transport rate (J_max) and the quantum efficiency for photon capture (α):

$$\begin{array}{}\text{(A4)}& J_{\mathrm{E}}=J\left(I\right){\displaystyle \frac{C_{\mathrm{i}}-{\mathrm{\Gamma }}_{\star }}{\mathrm{4}(C_{\mathrm{i}}+\mathrm{2}{\mathrm{\Gamma }}_{\star })}},J\left(I\right)=J_{\max }{\displaystyle \frac{\mathit{\alpha }I}{\sqrt{J_{\max }^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}{I}^{\mathrm{2}}}}}.\end{array}$$

A2 Soil water

In JSBACH the soil water budget is based on several reservoirs (e.g. skin, soil, bare soil, rain intercepted by canopy), and the different formulations are plentiful. We present here only the most crucial of these. Changes in volumetric soil water (θ_s, not to be confused with relative soil water content $\mathit{\theta }=\frac{{\mathit{\theta }}_{\mathrm{s}}}{{\mathit{\theta }}_{\mathrm{fc}}}$) due to rainfall (R), evapotranspiration (ET), snow melt (M), surface runoff (R_s) and drainage (D), are calculated with a geographically varying maximum field capacity (θ_fc) and soil water density (ρ_w).

$$\begin{array}{}\text{(A5)}& {\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{s}}}{\partial t}}=(\mathrm{1}-p_{\mathrm{int}})R+\mathrm{ET}+M-R_{\mathrm{s}}-D\end{array}$$

The interception parameter (p_int) also affects the amount of water intercepted by vegetation and bare soil, which further affects evaporation and transpiration. The skin reservoir is limited by w_skin, and excess water is transferred to soil water. Likewise when the soil water content (θ) is greater than parameter θ_dr, the excess water is rapidly drained (in addition to the limited drainage below this threshold), where d, d_min and d_max are constant parameters:

$$\begin{array}{}\text{(A6)}& D=d_{\min }\mathit{\theta }+(d_{\max }-d_{\min }){\left({\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{dr}}}{\mathrm{1}-{\mathit{\theta }}_{\mathrm{dr}}}}\right)}^d,\mathit{\theta }\ge {\mathit{\theta }}_{\mathrm{dr}}.\end{array}$$

Evaporation from wet surfaces (E_ws) depends on air density (ρ), specific humidity (q_a), saturation-specific humidity (q_s) at surface temperature (T_s) and pressure (p_s), and aerodynamic resistance (R_a). The aerodynamic resistance depends on heat transfer coefficient (C_h) and horizontal velocity (v_h).

$$\begin{array}{}\text{(A7)}& E_{\mathrm{ws}}=\mathit{\rho }{\displaystyle \frac{q_{\mathrm{a}}-q_{\mathrm{s}}(T_{\mathrm{s}},p_{\mathrm{s}})}{R_{\mathrm{a}}}},R_{\mathrm{a}}=C_{\mathrm{h}}\mathrm{|}{v}_{\mathrm{h}}{\mathrm{|}}^{-\mathrm{1}}\end{array}$$

Transpiration from vegetation (T_v) is likewise formulated but additionally depends on the stomatal resistance of the canopy (R_c), which is an inverse of the stomatal conductance, and as such, depends on which conductance model is used.

$$\begin{array}{}\text{(A8)}& T_{\mathrm{v}}=\mathit{\rho }{\displaystyle \frac{q_{\mathrm{a}}-q_{\mathrm{s}}(T_{\mathrm{s}},p_{\mathrm{s}})}{R_{\mathrm{a}}+R_{\mathrm{c}}}}\end{array}$$

Evaporation from dry bare soil (E_s) also has an added dependence on surface relative humidity (h_s) calculated from soil dryness:

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}& E_{\mathrm{s}}=\mathit{\rho }{\displaystyle \frac{q_{\mathrm{a}}-h_{\mathrm{s}}{q}_{\mathrm{s}}(T_{\mathrm{s}},p_{\mathrm{s}})}{R_{\mathrm{a}}}},\\ & h_{\mathrm{s}}=max\left[{\mathit{\theta }}_{\mathrm{hum}}(\mathrm{1}-\cos (\mathit{\pi }\mathit{\theta }\left)\right),min\left(\mathrm{1},{\displaystyle \frac{q_{\mathrm{a}}}{q_{\mathrm{s}}(T_{\mathrm{s}},p_{\mathrm{s}})}}\right)\right].\end{array}\end{array}$$

The total evapotranspiration is a weighted average of E_ws, T_v and E_s, where the weights are based on fill levels of reservoirs and the vegetative fraction of the grid cell.

A3 Logistic Growth Phenology model

The parameters from the LoGro-P are mainly used to determine the spring and autumn events for JSBACH. To determine the date of the spring event, we first introduce a few additional variables, namely the heatsum S_T(d), the number of chill days C(d) and the critical heatsum S_crit(d). T(d) denotes the mean temperature at day d.

$$\begin{array}{}\text{(A10)}& S_T\left(d\right)=\sum _{d^{\prime }=d_{\mathrm{0}}}^{d}max\left(T\right(d^{\prime })-T_{\mathrm{alt}},\mathrm{0})\end{array}$$

Heatsum S_T(d) cumulates the amount of “heat” above the parameter T_alt after the previous growing season. The actual starting date d₀ of the summation need not be known since it is enough to start the summation “reasonably late” after the last growth season.

$$\begin{array}{}\text{(A11)}& C\left(d\right)=\sum _{d^{\prime }=d_{\mathrm{a}}}^{d}H\left(T_{\mathrm{alt}}-T\left(d\right)\right)\end{array}$$

The number of chill days is calculated as the number of days when the mean temperature is below T_alt. Here H() denotes the Heaviside step function, and the summation starts at the day (d_a) of the last autumn event.

$$\begin{array}{}\text{(A12)}& S_{\mathrm{crit}}\left(d\right)=S_{\min }+S_{\mathrm{range}}{e}^{-C\left(d\right)/C_{\mathrm{decay}}}\end{array}$$

The critical heatsum (S_crit) decreases as the number of chill days C(d) increases, with an exponential memory loss parameter C_decay. The spring event happens when

$$\begin{array}{}\text{(A13)}& S_T\left(d\right)\ge S_{\mathrm{crit}}\left(d\right).\end{array}$$

The autumn event requires the definition of one more variable, the (pseudo) soil temperature (T_s(t)), which at time t is calculated as an average air temperature (T) with an exponential memory loss (T_ps). The autumn event occurs when T_s falls below a certain threshold. In the equation N is the normalisation constant and τ is the length of a time step.

$$\begin{array}{}\text{(A14)}& T_{\mathrm{s}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{N}}\sum _{n=-\mathrm{\infty }}^{t}T\left(n\right)e^{-(t-n)\frac{\mathit{\tau }}{T_{\mathrm{ps}}}}\end{array}$$

Appendix B: Stomatal conductance formulations

In this Appendix we present the stomatal conductance model formulations used in this study. In the original JSBACH formulation, the Baseline model (Knorr, 1997), the photosynthetic rate is resolved in two steps. First the stomatal conductance under conditions with no water stress is assumed to be controlled by photosynthetic activity (Schulze et al., 1994). Here the leaf-internal CO₂ concentration is assumed to be a constant fraction ($C_{\mathrm{i},\mathrm{pot}}=f_{\mathrm{C}\mathrm{3}}{C}_{\mathrm{a}}$) of ambient CO₂ concentration (C_a). This allows for an explicit resolution of the photosynthesis (Knorr, 1997). Then the impact of soil water availability is accounted for by a soil-moisture-dependent multiplier (β) that is identical for each canopy layer (Knorr, 1997).

$$\begin{array}{}\text{(B1)}& g_{\mathrm{s},\mathrm{pot}}={\displaystyle \frac{\mathrm{1.6}{A}_{\mathrm{n},\mathrm{pot}}}{C_{\mathrm{a}}-C_{\mathrm{i},\mathrm{pot}}}}\Rightarrow g_{\mathrm{s}}=\mathit{\beta }{g}_{\mathrm{s},\mathrm{pot}}\end{array}$$

After accounting for soil water stress, the net assimilation rate (A_n) and intercellular CO₂ concentration are (C_i) are recalculated using g_s and integrated over the leaf area index to produce canopy level estimates.

In the Bethy approach (Knorr, 2000), the unstressed canopy conductance (G_c,pot) is calculated similarly to the Baseline model but is potentially further limited by the water supply function of the maximum transpiration rate (T_supply=βT_max). T_max is a fixed and predefined upper limit for transpiration as in Knauer et al. (2015).

$$\begin{array}{}\text{(B2)}& \begin{array}{rl}& G_{\mathrm{c}}=\left\{\begin{array}{ll}{G}_{\mathrm{c},\mathrm{pot}}\frac{T_{\mathrm{supply}}}{T_{\mathrm{pot}}},& T_{\mathrm{pot}}\ge T_{\mathrm{supply}}\ge \mathrm{0}\\ G_{\mathrm{c},\mathrm{pot}},& T_{\mathrm{pot}}<T_{\mathrm{supply}}\end{array}\right.\\ & T_{\mathrm{pot}}=\mathit{\rho }{\displaystyle \frac{q_{\mathrm{s}}-q_{\mathrm{a}}}{\mathrm{1}/G_{\mathrm{a}}+\mathrm{1}/G_{c,\mathrm{pot}}}}\end{array}\end{array}$$

The potential (unstressed) transpiration rate (T_pot) is a function of air density (ρ), saturation-specific humidity (q_s) at given temperature and pressure, specific humidity (q_a), aerodynamic conductance (G_a), and unstressed canopy conductance (G_c,pot). After this scaling, the net assimilation rate and intercellular CO₂ concentration are recalculated as in the Baseline model.

The Ball–Berry variants relate the stomatal conductance (g_s) to empirically fitted parameters g₀ (mol m⁻² s⁻¹) and g₁ (unitless, except for $g_{\mathrm{1}}^{\mathrm{USO}}$, which has units of $\sqrt{\text{kPa}}$) that respectively represent the residual stomatal conductance and the slope of the function. The stomatal conductance is a function of the net assimilation rate (A_n), the water stress factor (β) and the atmospheric CO₂ concentration (C_a). The original Ball–Berry formulation (Ball et al., 1987) also depends on relative humidity at leaf surface (h_s). In the Leuning model (Leuning, 1995), the CO₂ concentration is reduced by the CO₂ compensation point (Γ) as well as scaled by the vapour pressure deficit (D_s) and a constant (D₀) depicting the stomatal sensitivity to changes in D_s. The Friend and Kiang model (Friend and Kiang, 2005) adds an exponential dependency on the difference of specific (q_a) and saturation-specific humidity (q_sat) with empirically fitted constants a=2.8 and b=80. The unified stomatal optimisation model (Medlyn et al., 2011) also adds a dependency to the vapour pressure deficit (D_s).

$$\begin{array}{}\text{(B3)}& {\displaystyle }\begin{array}{rl}& g_s^{\mathrm{BB}}=g_{\mathrm{0}}^{\mathrm{BB}}+g_{\mathrm{1}}^{\mathrm{BB}}\mathit{\beta }{\displaystyle \frac{A_{\mathrm{n}}{h}_{\mathrm{s}}}{C_{\mathrm{a}}}}\\ & g_{\mathrm{s}}^{\mathrm{Leu}}=g_{\mathrm{0}}^{\mathrm{Leu}}+g_{\mathrm{1}}^{\mathrm{Leu}}\mathit{\beta }{\displaystyle \frac{A_{\mathrm{n}}}{(C_{\mathrm{a}}-\mathrm{\Gamma })(\mathrm{1}+D_{\mathrm{s}}/D_{\mathrm{0}})}}\\ & g_{\mathrm{s}}^{F\mathit{\&}K}=g_{\mathrm{0}}^{F\mathit{\&}K}+g_{\mathrm{1}}^{F\mathit{\&}K}\mathit{\beta }{\displaystyle \frac{A_{\mathrm{n}}{a}^{-d(q_{\mathrm{sat}}-q_{\mathrm{a}})}}{C_{\mathrm{a}}}}\\ & g_{\mathrm{s}}^{\mathrm{USO}}=g_{\mathrm{0}}^{\mathrm{USO}}+\mathrm{1.6}\left(\mathrm{1}+{\displaystyle \frac{g_{\mathrm{1}}^{\mathrm{USO}}\mathit{\beta }}{\sqrt{D_{\mathrm{s}}}}}\right){\displaystyle \frac{A_{\mathrm{n}}}{C_{\mathrm{a}}}}\end{array}\end{array}$$