Constraining ecosystem model with Adaptive Metropolis algorithm using boreal forest site eddy covariance measurements

We examined parameter optimization in JSBACH ecosystem model, applied for two boreal forest sites (Hyytiälä and Sodankylä) in Finland. We identified and tested key parameters in soil hydrology and forest water and carbon exchange related formulations and optimized them using the Adaptive Metropolis algorithm (AM) for Hyytiälä with a five year calibration period (2000–2004) followed by a four year validation period (2005–2008). Sodankylä acted as an independent validation site, where optimizations were not made. 5 The tuning provided estimates for full distribution of possible parameter, along with information about correlation, sensitivity and identifiability. Some parameters were correlated with each other due to phenomenological connection between carbon uptake and water stress or other connections due to the set-up of the model formulations. The latter holds especially for vegetation phenology parameters. The least identifiable parameters include phenology parameters, parameters connecting relative humidity and soil dryness, and the field capacity of the skin reservoir. These soil parameters were masked by the large 10 contribution from vegetation transpiration. In addition to leaf area index and maximum carboxylation rate, the most effective parameters adjusting GPP and ET fluxes in seasonal tuning were related to soil wilting point, drainage and moisture stress imposed on vegetation. For daily and half-hourly tunings the most important parameters were the ratio of leaf internal CO2 concentration to external CO2 and the parameter connecting relative humidity and soil dryness. Effectively the seasonal tuning transferred water from soil moisture into ET, and 15 daily and half-hourly tunings reversed this process. The seasonal tuning improved the month-to-month development of GPP and ET, and produced the most stable estimates of water use efficiency. When compared to the seasonal tuning, the daily tuning is worse on the seasonal scale. However, daily parametrization reproduced the observations for average diurnal cycle best, except the GPP for Sodankylä validation period, where half-hourly tuned parameters were better. In general, the daily tuning providing the most reduction in model20 data mismatch. The models response to drought was unaffected by our parametrizations and further studies are needed into enhancing the dry response in JSBACH.


Introduction
Inverse modelling of ecosystem model parameters against in situ observations is an established way to tune model parameters (e.g.Scharnagl et al., 2011).As observation sites have their own characteristics, it is necessary to make local site simulations for model evaluation and calibration as they may reveal new insight into model behaviour and guide further development.Model-data fusion has been applied for boreal forest sites by, e.g., Aalto et al. (2004) Peltoniemi et al. (2015b), Thum et al. (2007Thum et al. ( , 2008) ) and Wu et al. (2011).
Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
In this study we perform site level parameter optimisation in the JSBACH model (Kaminski et al., 2013;Knorr and Kattge, 2005;Reick et al., 2013).JSBACH is the land surface component of the Earth system model of the Max Planck Institute for Meteorology (MPI-ESM), used to simulate water and carbon storages and fluxes in the soil-vegetationatmosphere continuum.The water and carbon fluxes are coupled in the model and thus modification of parameters related to one component usually has an effect on the others as well.The optimisation process and the optimised values are also affected by the assimilation frequency and interval in minimising the model-data mismatch.This effect can be studied in numerous ways; e.g.Santaren et al. (2014) varied the length of assimilation interval while Matheny et al. (2014) focused on the diurnal error patterns.
The motivation for this study comes from results showing that CMIP5 model simulations, one of which is MPI-ESM, have systematic evapotranspiration biases over continental areas (Mueller and Seneviratne, 2014).These kinds of biases not only have significant implications for climate change projections (Boé and Terray, 2008) but also have a distinctive behaviour on a regional scale.In addition, a comparative study of the gross primary production (GPP) of Finnish forests (Peltoniemi et al., 2015a) revealed that JSBACH has an insufficient response to water limitation in Finland -it tends to overestimate GPP and evapotranspiration during dry periods.This is especially apparent in the dry year 2006, as JSBACH is unable to transfer the reduced rainfall into lower levels of GPP.
In this study we apply the JSBACH ecosystem model for Hyytiälä (Kolari et al., 2009;Suni et al., 2003) and Sodankylä (Aurela, 2005;Thum et al., 2008) sites.We identify key parameters in soil hydrology and evapotranspirationrelated formulations and test their effectiveness with elementary methods.We study the effect of different timescales (seasonal, daily and half-hourly) on the assimilation process and the effect of this on the optimised parameter values.Several optimisations are performed using the adaptive Metropolis (AM) algorithm over a 5-year calibration period (2000)(2001)(2002)(2003)(2004) followed by a 4-year validation period (2005)(2006)(2007)(2008).
The goals of this study are to test the applicability of the AM optimisation method for JSBACH and the impact of different temporal resolutions on the optimisation process, and to improve the models response to environmental drivers, focusing on dryness.

Measurements, sites and instrumentation
In this study we use site level data from two Finnish measurement sites: Hyytiälä (61 • 51 N, 24 • 17 E; 180 m a.s.l.) and Sodankylä (67 • 22 N, 26 • 38 E; 179 m a.s.l.).These wellestablished sites have long continuous measurement data sets representing the southern and northern boreal Finnish coniferous evergreen forests.The data used in this study are available for the scientific community through various databases such as FLUXNET (re3data.org, 2016).
Hyytiälä site is a Finnish Scots pine (Pinus sylvestris) forest (Kolari et al., 2009), planted in 1962 after burning and mechanical soil preparation.The soil type in Hyytiälä is Haplic Podzol on glacial till and the site is of medium fertility (Kolari et al., 2009).The forest also has sparse understory of Norway Spruce (Picea abies) and scattered deciduous trees.The maximum of measured all-sided leaf area index (LAI) is 6.5 m 2 m −2 for the Scots pine.The carbon dioxide (CO 2 ) and water vapour (H 2 O) fluxes between vegetation and atmosphere have been measured in Hyytiälä continuously since 1997 (Vesala et al., 2005).
The Sodankylä forest, in Sodankylä at the Finnish Meteorological Institute's Arctic Research Centre, is also a Scots pine forest (Pinus sylvestris) with maximum measured LAI of 3.6 m 2 m −2 as determined from a forest inventory in early autumn (Thum et al., 2007).The forest on fluvial sandy Podzol has been naturally regenerated after forest fires with tree age ranging approximately from 50 to 100 years.The sparse ground vegetation consists of lichens (73 %), mosses (12 %) and ericaceous shrubs (15 %).The CO 2 and H 2 O flux measurements in Sodankylä have been running since 2000 (Aurela, 2005).
The CO 2 and H 2 O fluxes were measured by the micrometeorological eddy covariance (EC) method, which provides a direct measurement of the mass and energy exchange between the atmosphere and the biosphere averaged on an ecosystem scale.In the EC method, the flux is obtained as the covariance of the high-frequency (10 Hz) observations of vertical wind speed and the constituent in question (Baldocchi, 2003).The CO 2 fluxes were corrected for the storage change below the measurement height to accurately estimate the net ecosystem CO 2 exchange (NEE).The GPP was derived by subtracting the modelled respiration (R) from the NEE observation (GPP = NEE − R) utilising standard flux partitioning procedures (Reichstein et al., 2005;Kolari et al., 2009).By using the same parametrisations as in the partitioning, the NEE and GPP time series were gap-filled for comparison with the model results.The daily evapotranspiration (ET) sums were calculated from H 2 O flux data that were gap-filled based on the mean diurnal cycles or regressions on available radiative energy.
The EC instrumentation in Hyytiälä consisted of a Solent 1012R3 three-axis sonic anemometer (Gill Instruments Ltd., Lymington, UK) and a LI-6262 closed-path CO 2 /H 2 O gas analyser (Li-Cor Inc., Lincoln, NE, USA), while in Sodankylä a USA-1 (METEK GmbH, Elmshorn, Germany) anemometer and an LI-7000 (Li-Cor., Inc., Lincoln, NE, USA) closed-path gas analyser were used.The EC fluxes were calculated as half-hourly averages taking into account the required corrections.The measurement systems and the post-processing procedures have been presented in more detail for Hyytiälä by Kolari et al. (2004) and Mammarella et al. (2009), and for Sodankylä by Aurela (2005) and Aurela et al. (2009).
The measurement error in the EC flux data may be classified into two categories: systematic errors and random errors.The main systematic errors (density fluctuations, highfrequency losses, calibration issues) are mostly corrected for 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 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).

The JSBACH model
JSBACH is a process-based ecosystem model and the land surface component of the MPI-ESM.We used JSBACH offline using an observational atmospheric data set 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 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 (a single grid point), the grid resolution does not affect the results (Tesfa et al., 2014;Singh et al., 2015).We give here a general introduction to JSBACH, whereas a more complete model description can be found in Roeckner et al. (2003).
In JSBACH the land surface is a fractional structure where the land grid cells are divided into tiles representing the most prevalent vegetation classes called plant functional types (PFTs) within each grid cell (Reick et al., 2013).The grid cell is first divided into bare soil and vegetative area which is furthermore fractionally divided into PFTs.The model was set up to effectively use only one tile, coniferous evergreen trees.Henceforth, all model and process descriptions are considered in relation to coniferous evergreen trees and no distinction between PFT-specific and general parameters are made in this study.
Coniferous evergreen trees are characterised by a set of parameters that control vegetation-related biological and physical processes accounting for the land-atmosphere interactions.We made use of expert knowledge to set these parameters for our local sites and verified that they are in line with those presented by Hagemann (2002) and Hagemann and Stacke (2015).
The seasonal development of LAI is regulated by air temperature and soil moisture with a specific maximum LAI as a limiting value.The cycle is driven by a pseudo soil tem-perature that is a weighted running mean of air temperature.The predictions of phenology are produced by the Logistic Growth Phenology (LoGro-P) model of JSBACH.
Photosynthesis is described by the biochemical photosynthesis model (Farquhar et al., 1980).Following Kattge et al. (2009), we set the maximum carboxylation rate at 25 • C to 1.9 times the maximum electron transport rate at 25 • C.
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 2 concentration is assumed to be a constant fraction of ambient concentration, which 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).
Radiation absorption is estimated by a two-stream approximation within a three-layer canopy (Sellers, 1985).Especially in the sparse canopies, the radiation absorption is affected by clumping of the leaves, which is here taken into account according to the formulation by Knorr (1997).

The JSBACH model spin-up and runs
Before tuning the JSBACH model, some of the more slowly changing variables (e.g.LAI) need to be equilibrated in order to bring the model into a (semi-)steady state.We achieve this by running the model through a spin-up period generated by looping the measurement interval over itself.During this period the necessary variables are equilibrated and their values become acceptable for the tuning process.At the end of the spin-up a restart file is generated so that the model can be restarted from this state.
We use half-hourly measurements from 1999 to 2008 for Hyytiälä.The spin-up finishes at the end of 1999 and is followed by a calibration period (abbreviated as HC for Hyytiälä calibration) of 2000-2004 and a validation period (HV for Hyytiälä validation) of 2005-2008, including an exceptionally dry summer in 2006.The set-up for Sodankylä is similar but we use measurements from 2000 to 2008, where the spin-up finishes at the end of 2008.The model is then restarted from the start of 2000, but we only examine the Sodankylä validation period (SV) of 2005-2008.The main reason to exclude the Sodankylä calibration period is that essentially we do not calibrate (or tune) the model for Sodankylä and we do not want to appear to do so.
The meteorological data used to drive the climate were air temperature, air pressure, atmospheric CO 2 concentration, precipitation, specific humidity, short-and longwave radiation, potential shortwave radiation and wind speed.

The parameters
The JSBACH model was modified to fit our custom-built test bed so that all parameters of interest could be read from an Eq.(A6) * These parameters were tested but yielded no or only a minimal response to cost functions and were thus removed from the trial.
external file.We examined 15 parameters (Table 1) that are for convenience separated into three classes.The class I parameters are used differently from those of class II and IIInamely, class I parameters are only tuned in the seasonal tuning (explained in detail in Sect.3.1).Additionally, the only distinction between class II and III parameters is that the latter belong to a specific part of JSBACH called the LoGro-P -there is no difference in how these parameters are used.We also note that the only parameter (of those examined) that can vary from site to site is veg max (the vegetative fraction of a grid cell).

Parameter sampling
The parameter sampling in this study was done with the AM algorithm.The AM algorithm is an adaptive Markov chain Monte Carlo (MCMC) process described below (it is not strictly Markovian but satisfies the necessary ergodicity requirements).AM is based on the classical Metropolis algorithm, extended with the adaptation of the parameter proposal distribution.Due to the adaptive nature of AM, it does not rely on the choice of the initial proposal distribution.AM is a sampling method that produces estimates of the full distribution of possible parameter values (unlike straightforward optimisation methods), thus enabling, e.g., the study of parameter identifiability, sensitivity and (nonlinear) correlation -this information is paramount to understanding the optimisation process in contrast to merely receiving the optimised parameter values.The rigorous mathematical presentation of the AM algorithm is given in Haario et al. (2001).
The AM algorithm draws samples (sets of parameters) from the parameter space to generate probability distributions for the parameters.The consecutive draws form an MCMC chain.We used the algorithm simultaneously for several independent chains that are parallel adaptations of the algorithmic process (see e.g.Craiu et al., 2009;Solonen et al., 2012) -we take several random starting points and launch the algorithm for each of these simultaneously.The history of all chains is used for updating the proposal covariance matrix that describes how the parameters relate to one another.Our initial proposal covariance matrix had diagonal elements corresponding to 1/200 of the respective parameter's range.The first sample for each chain was chosen at random within this range.The algorithmic process can be described by a few steps: 1. Draw a new sample (x ) of the parameter space from the vicinity of the current sample (x) using the current proposal covariance matrix.
2. Calculate the acceptance ratio (a) for the drawn sample.This is the value of a likelihood function (f ) that is proportional to the desired probability distribution, at the drawn sample divided by the value at the current sample (a = f (x )/f (x)).
3. Accept the new candidate (x ) with the probability a (if a ≥ 1, we always accept).If the candidate was rejected, the current sample (x) is reused as a basis of the next draw and repeated in the chain.Update the covariance matrix and draw a new sample.
We obtain the likelihood function (f ) from the cost functions (cf) described below by assuming Gaussian error statistics and setting f = e −cf .In general to estimate the distribution of parameters of any model based on some data, we require some information about the underlying measurement and modelling errors.We treat the JSBACH model as described by the equation y = M(x, θ ) + e.Here y are the observations, x is the model state vector, θ are the current parameters and e is the model-data mismatch.Since we only have a robust estimate for the measurement errors and no true error statistics for the model, the full error (e) is treated as Gaussian white noise.
The cost function (Eq. 1) used in this study for seasonal tuning is based on summary statistics of GPP and evapotranspiration (ET) along with the maximum of LAI.The cost function (Eq. 1) calculates the relative error in total GPP, ET and growing season maximum of LAI against observations (these are respectively denoted as G 1 , E 1 and L 1 ) and sums them up.Overlined variables refer to the mean value of that variable for a given period (calibration or validation), subscripts denote observation or model results.
The second cost function (Eq.2) is a slightly modified mean squared error estimate used for daily (cf 2 ) and halfhourly (cf 3 ) tuning.With multiple variables there is always the problem of having one variable dominating over the others.Since no true errors were available, it was decided to normalise the residuals using the mean of observations in the cost function (Eq.2).This way the resulting function is sensitive to changes in both variables -AM is used as a noiseresistant optimiser and sampling is done in the spirit of studying the identifiability and correlations of the parameters.The components are denoted as G 2 , E 2 for daily and G 3 , E 3 for half-hourly tuning.
As noted previously, JSBACH was used uncoupled from the other components of the full MPI-ESM.This has a tendency to lead to biased results in the model runs as has been recently studied by Dalmonech et al. (2015).Especially in the high latitudes, evapotranspiration can be unrealistic during winter since night-time is longer and temperatures low.In order to improve the credibility of our results, we masked the evapotranspiration values of the coldest periods, and only took into account those from May to September for each year in both cost functions.

Parameter analysis
The optimised parameter values are taken as the mean values of all chains in the sampling process.In the case that the parameter chains converge to a bound of an a priori prescribed range of allowed values, the maximum a posteriori (MAP) value is used instead.After tuning the model, we analysed different aspects of this process.Class I parameters are excluded from this analysis since they are used to bring the model to an "acceptable initial state"; hence, we regard them as a part of the model initialisation (our motivation is explained in Sect.3.1).We calculated the correlations and correlation matrices between different parameters for different tunings using the tested parameter vectors in the AM process.Then we performed a principal component analysis (PCA) on the correlation matrices to get the eigenvectors (v i ) and eigenvalues (e i ) of the least identifiable parameters in the tuning process with the given data.The PCA transforms the correlation matrix into an orthogonal form where the eigenvector related to the greatest eigenvalue is the least identifiable with the given data.We then calculate the weight (w i = e 2 i i e 2 i ) for each component (or vector v i ; note that the squared weights sum up to one).We also determine the most dominant parameters for each component (v i ) by similarly dividing the length of the vector towards that parameter by the length of the whole vector (weight of vector components).
The information derived with PCA could be extracted by analysing the parameters posterior probability distributions, but PCA yields a simple, straightforward method for the same purpose.The main caveat of the standard PCA method is that it is not applicable to cases with strong nonlinear correlations.Therefore, we also calculate kernel density estimates (KDE) for the parameters to show that there are no nonlinear correlations.The KDE method places a Gaussian distribution (kernels) centred at each parameter of the MCMC chain and then sums these kernels to produce an estimate for the whole distribution.The bandwidth is calculated using the Scott's rule (Scott, 2004).
We also wanted to examine which parameters contributed the most to the change in the cost function values as we switched from one parameter set to another.This was done by calculating the change in the cost function values of the tuned parameter set and a set where one parameter has been reverted to the value the tuning started with (henceforth, the reference values are for seasonal tuning the default values www.nonlin-processes-geophys.net/23/447/2016/ Nonlin.Processes Geophys., 23, 447-465, 2016 and for daily and half-hourly tuning the seasonally tuned values).We call this method "relative effectiveness", since we want to analyse the effect of the parameters to the cost function.For each tuned set of parameter values, the relative effectiveness of a parameter is calculated as follows: 1. change one parameter from the set of tuned parameter values to a reference value and calculate the difference in the cost function for the changed set and the tuned set; 2. return the changed parameter to the tuned value and repeat for all parameters (sum up the differences); 3. the relative effectiveness for each parameter is the difference obtained from step 1 divided by the sum from step 2.
The relative effectiveness is similar to a class of methods commonly referred to as the one-at-a-time (OAT) or onefactor-at-a-time (OFAT) methods.These methods are generally used to acquire robust information about model behaviour when one parameter at a time is changed to a new and hopefully better value (e.g.Murphy et al., 2004).The main difference of our method to classical methods such as the Morris OAT (Morris, 1991) is that in such methods the change in values is (usually) random, whereas we have fixed values.Additionally, our point of view is from the optimised parameters to the original state -we have already optimised the parameters (as a group) and merely want some robust and easily comprehensible information about the effect of changes in parameter values to the cost functions.This method does not reveal information about how well the parameters constrain the cost function (e.g.we could have a highly dominating parameter that would optimise to the default value and hence the relative effectiveness would be zero), rather which parameters contribute most to the change in cost function values.Lastly, we calculate the root mean squared error (RMSE; 2 ) for the time series generated by the different tunings (o i is observed and m i is modelled).

Model tuning
The model was optimised for Hyytiälä with the AM algorithm using three different timescales: seasonal, daily and half-hourly tuning, which are described below.Tuning was done on a powerful laptop with an Intel Core i7-3520M processor.We removed unwanted output streams from the model and tweaked the I/O.With a single core the spin-up generation takes approximately 150 s, the run through calibration period with daily output takes 20 s and with half-hourly output 320 s.We used daily output also for the seasonal tuning.

Seasonal tuning
The fundamental motivation for the seasonal tuning is to ensure that the model reproduces the observed growing season maximum of LAI, since we have previously noticed that JSBACH underestimates LAI at the site level (even the default value of max is lower than the measured maximum for Hyytiälä).The reason for this approach was to enhance the vegetation transpiration and to emphasise the model response to precipitation.We also want to ensure that the model performs adequately well in terms of seasonal cumulative GPP and ET.The seasonal tuning was done in three consecutive steps each using the cost function (Eq.1).The procedure is as follows: 1.All three class I parameters are tuned with four independent chains each consisting of 3000 samples.This step required a 30-year spin-up for each sample separately.
2. Class II and III parameters are each separately tested with 24 evenly separated values for an extensive range and those nine parameters that did not yield a negligible difference in the maximal and minimal values in the objective function are tuned.The consequent tuning was done with eight independent chains each consisting of 10 000 samples.A single spin-up, common for all samples, used optimal parameter values from step 1 and default values for the rest of the parameters.
3. All the previously tuned 12 parameters with eight independent chains each consisting of 10 000 samples are returned.Initial proposal covariance was generated from previous step and spin-up was generated separately for each sample.
At the end of seasonal tuning, class I parameters were fixed and a single spin-up was generated to be used with daily and half-hourly tuning.This approach is computationally justifiable (as we do not have to rerun the spin-up at each iteration of the algorithm) and is also acceptable from a modelling point of view since the robust site level scaling has already been done.The vegetative fraction of a grid cell remained at its default value of 0.52 and the carboxylation rate at 25 • C was lowered to 45.0 (and the electron transport rate to 85.5).

Daily and half-hourly tuning
The difference in daily and half-hourly tuning is the time interval used in the model output and observations in the cost function (Eq.2).For both tuning runs we first tested the response of class II and III parameters against the cost function (Eq.2) and removed those parameters that yielded only negligible or no response (as in step 2 in "Seasonal tuning").The rest of the parameters (12) were then tuned using eight independent chains each consisting of 10 000 samples.
It should be noted that even though the cost function (Eq.2) formulation is the same for daily and half-hourly tun- ing, the values of the cost function are not directly comparable.Half-hourly tuning uses 48 values per day, and the resulting diurnal pattern resembles the form of the normal distribution.In daily tuning we use an average of these values.In practice, the component and cost function values will be higher for half-hourly tuning.

Tuning for Sodankylä
After tuning the model for Hyytiälä we took the parameter set from seasonal tuning and re-tuned only the maximum LAI parameter ( max ) with the cost function (Eq. 1) for Sodankylä.This was done because the measured LAI for Sodankylä is approximately half of that for Hyytiälä.The other parameter values were taken from the respective Hyytiälä tuning runs and spin-ups were generated similarly to Hyytiälä spin-ups so that we could use the Sodankylä results to validate the tuning process.

Results and discussion
The parameters and cost function components involved in the different phases of the optimisation process need to be studied before the performance of the optimisation method can be evaluated.
As noted above, we decided to reject the unreliable wintertime ET values.This masking leaves out the start of the growing season, which reduces the reliability of some of the tuned parameters, including all the LoGro phenology model parameters (class III), which mostly affect the timing of the spring event and regulate the development of the LAI towards the peak season.However, as a result of the tuning processes, all the analysed parameters were revealed to have unimodal posterior probability distributions, with different skewness and deviations.
We analysed the correlations and effectiveness of the parameters in the seasonal, daily and half-hourly optimisations on the Hyytiälä site for the calibration period.We also analysed the contributions from the cost function components referring to ET, GPP and LAI and generated the time series Table 3. Significant components of principal component analysis for the different tunings.The given parameters are the most dominant within the component and the ratio is how many times larger the factor related to the first parameter is when compared to that of the second.Coverage reveals how much of the component is accounted for by the given parameters (sum of the weights of given vector components).and daily cycles of GPP and ET for both Hyytiälä and Sodankylä sites.For all these examinations, individual spin-ups were generated using the optimised parameter values.

Component
The parameter correlations (Table 2) do not reveal much information, which is common for larger systems where the underlying parameter dependencies are more complex.Usually more sophisticated methods are used to analyse the parameters, but we omit these examinations here since pairwise Kernel density estimates (Fig. 1) did not reveal any new insights.
The strongest correlation was between the ratio of leaf internal CO 2 concentration to external CO 2 (f C3 ) and fraction of soil moisture above which transpiration is unaffected by soil moisture stress (w tsp ) in all the tunings.This positive correlation strengthens as we increase the temporal resolution (and the complexity of the underlying cost function).This is due to the carbon assimilation that is limited not only by the amount of carbon available but also by a linear water stress factor (which takes the value of zero at the wilting point (w wilt ) and one at the w tsp ), which is checked at each time step.Most of the other parameters with high correlations are those of the LoGro phenology model, where we would expect high correlation since the parameters are intimately connected.
Approximately half of the parameters with high correlation are also the least identifiable (Table 3) with the given data and cost function.This means that the values these parameters acquire, as a result of the tuning process, are the most unreliable -it does not reflect on the parameters contribution to the cost function.The PCA merely highlights where most of the parametric unreliability lies.
The PCA analysis revealed that most of the unreliability is explained by a handful of parameters.Disregarding those of the LoGro phenology model, the two most dominantly unreliable parameters in every tuning were the fraction depicting relative humidity based on soil dryness (w hum ) and the maximum field capacity of the skin reservoir (w skin ).Both of these parameters affect the amount of water available for evaporation from bare soil and are both subject to changes in other parameters.Bare soil evaporation is also dominated by vegetative transpiration, which explains why these two parameters are the most unreliable.

The parameters and their relative effectiveness
The default and optimised parameter values from the different tuning metrics are presented in Table 4 along with their relative effectiveness.The reference values for seasonal tuning are the default values.Since we fixed class I parameters with seasonal tuning, the realistic reference values for daily and half-hourly tunings are the seasonal parameter val- ues.Here we note that using one spin-up for all daily and half-hourly optimisation runs is computationally justifiable but generates errors as the general spin-up differs from those generated by the optimised parameters.These errors are relatively small but give rise to, e.g., the negative relative effectiveness values in daily and half-hourly parametrisations.
Most seasonally tuned parameters are near their default values and the most effective parameters are the fraction of soil moisture above which transpiration is unaffected by soil moisture stress (w tsp ), the fraction of soil moisture at permanent wilting point (w pwp ) and the fraction of field capacity above which fast drainage occurs (w dr ).For daily and halfhourly tunings the most important parameters are the ratio of leaf internal CO 2 concentration to external CO 2 (f C3 ) and the fraction depicting relative humidity (w hum ).It should be noted that w hum was one of the least identifiable parameters for seasonal tuning.Taking into account the importance of these parameters on transpiration and soil moisture estimations, we took a closer look at modelled soil moisture and evapotranspiration components for the calibration period (taking into account only values from May to September for each year as explained at the end of Sect.2.5. When we compare the model output streams with seasonal against those with default parametrisation, we notice that the average evapotranspiration for the calibration period has increased 15 %.Most of this is due to not only added transpiration (18 % increase) but also increased evaporation (6 %).In addition drainage was accelerated by 11 %.These increases are mostly compensated by a 15 % reduction in average soil moisture.In addition soil moisture values that are under the limit when transpiration is affected by soil moisture stress (below the value of w tsp ) increased 2.3 %.
The daily and half-hourly tunings lower the average evapotranspiration by 22 and 35 % respectively, when compared to the seasonal values.Transpiration is decreased by 28 and 37 %, whereas evaporation is increased by 0.5 % and de-creased by 28 % respectively, for daily tuning and half-hourly tuning.Soil moisture is increased by 11 and 8 % and the amount of values below w tsp is decreased by 62 % for daily tuning and increased by 7 % for half-hourly tuning.Out of curiosity, both the adjustment parameter in stability functions (c b ) and the fraction of precipitation intercepted by canopy (p int ) have been significantly increased with daily tuning and returned to seasonally tuned values with halfhourly tuning.

The cost function components
Using the optimised values (parametrisations), we calculated the components of each cost function for Hyytiälä calibration period and Hyytiälä and Sodankylä validation period (Table 5).
First, we note that with the default parameters L 1 dominates cf 1 for Hyytiälä and contributes approximately 90 % to its value.As expected the L 1 for Sodankylä is not as dominant as for Hyytiälä since the measured maximum of LAI for Hyytiälä is roughly half as large as for Sodankylä, which directly lowers the LAI component in cost function (Eq.1).The L 1 contribution is significantly reduced with the seasonally tuned parameters as was our intention and even though LAI plays no part in daily and half-hourly tunings, the differences in the maximum value are negligible.
Second, the value of the E 1 component (error in seasonal ET) with default parametrisation is significantly increased in daily and especially half-hourly parametrisations.Simultaneously the value of G 1 is significantly lowered.The component values for seasonal parametrisation are better than the default values with the exception of E 1 for Hyytiälä validation period.
Third, for the cost function (Eq.2) the pairwise ratio of dominating E i or G i components in all tunings is 5 : 1.On average E 2 /E 3 contributes to approximately 60% of cf 2 /cf 3 .
www.nonlin-processes-geophys.net/23/447/2016/ Nonlin.Processes Geophys., 23, 447-465, 2016 Table 5. Cost function components for each parametrisation for Hyytiälä calibration (HC), Hyytiälä validation (HV) and Sodankylä validation (SV) periods.L 1 , E 1 and G 1 are the LAI, ET and GPP components in cost function (Eq.1), represented by cf 1 and used for seasonal tuning.Likewise E 2 and G 2 are the components in cost function (Eq.2) for daily values (cf 2 ), whereas E 3 and G 3 are for half-hourly values (cf 3 ).Note that the values of cf 2 and cf 3 are not directly comparable.This translates to ET being twice as significant as GPP in the cost function (Eq.2).The main reason for ET dominating GPP is that ET is more erratic in comparison to GPP and the residuals of ET (divided by the mean value) are larger than the residuals of GPP.The daily and half-hourly tunings themselves work as intended as they lower the corresponding cost function value.It is noteworthy to mention that the G 2 component gets its lowest value for both validation periods with the half-hourly parametrisation even though G 2 calculates GPP errors on a daily scale.Lastly, we examine how the algorithm and cost functions have performed.The best parameter set (the lowest cost function value) for a given cost function, in each of the three different periods (HC, HV, SV), is the same as that used in the corresponding tuning process.For example the lowest value for cf 1 (the cost function for seasonal tuning) in Sodankylä validation period (0.07) coincides with the seasonally tuned parameters.This is expected as the tuning process aims to be the "best" parameter value, which reassures us that no gross mistakes (human errors) have been made.The relation holds true for every cost function with the exception of cf 1 for Hyytiälä validation period, where the lowest value is reached with the daily tuned parameters (we note that the absolute difference between daily and seasonally tuned parameters is small).Hence we can confidently state that the algorithm and cost functions have performed as intended, especially since the optimised parameters work for Sodankylä as well, where no optimisation (besides the site-specific maximum of LAI) was applied.

Time series
The overall structure of the model time series was not affected by the parametrisations obtained with different tunings (Figs. 2 and 3).Some time series characteristics have been enhanced and others reduced but the timing of the peaks and dips in GPP and ET are the same as before.The corresponding RMSE and bias estimates are given in Table 6.In comparison we estimated the PRELES model biases for Hyytiälä from Fig. 5 in Peltoniemi et al. (2015b).These estimates give a bias of 0.81 × 10 −6 kg m −2 s −1 (0.07 mm m −2 day −1 ) for ET and −1.45 × 10 −7 mol [CO 2 ] m −2 s −1 (−0.15 g(C) m −2 day −1 ) for GPP.Additionally, the coefficient of determination (r 2 ) for GPP in Hyytiälä is in the range of 0.74-0.76for all tunings, whereas the values reported in literature range from 0.68 (Trusilova et al., 2004) to 0.96 (Peltoniemi et al., 2015b) with most above 0.9 (Aalto et al., 2004;Duursma et al., 2009).For additional comparisons see also Abramowitz et al. (2007).Note that our estimates are calculated using only values from the beginning of May to the end of September.
The best seasonal performance was obtained by seasonal tuning as we previously noted from the cost function components (Table 5).Even though the optimisation is done on the seasonal level, especially the GPP cycle is noticeably improved from that generated by the default parameters.This tuning also gives rise to the most stable (least fluctuating) water use efficiency (WUE), when calculated as a pointwise ratio of GPP and ET.We use WUE here only as a diagnostic variable to examine the balance between the GPP and ET.
When compared to the seasonal tuning, the daily tuning is worse on the seasonal scale and lowers both the ET and GPP cycles.WUE follows the observations better but starts to give rise to some fluctuation.With half-hourly tuning, this behaviour is further enhanced and especially ET is lowered to too low levels, which manifests the high WUE values.The Table 6.RMSE and bias of ET and GPP calculated from half-hourly data for first two summers of the validation period for Hyytiälä (corresponding to Fig. 2) and last two summers of the validation period for Sodankylä (corresponding to Fig. 3).worsening in the model time series with daily and half-hourly tunings are explained by biases in the diurnal cycle.

Diurnal cycles
Average diurnal cycles with different parametrisations (Fig. 4) show that modelled night-time ET values are too low when compared to the observed and this behaviour was not affected by the tunings.Low night-time values are compensated by too high midday values in the default and seasonal tuning so that the average daily and seasonal values are on an acceptable level.For the daily and half-hourly tuning, the algorithm lowers the daytime values, which results in too low average daily and half-hourly values.It is noteworthy to mention that with the default setting we get too low GPP for Hyytiälä but too high GPP for Sodankylä.The unrealistic wintertime and the biased night-time ET values actually have the same origin.Since we do not have the coupling from the land surface model (LSM) back to the atmosphere, we get an erroneous energy balance as we lose the energy released by condensation.Disregarding the default parametrisation we notice that seasonal parametrisation show the highest values, daily in the middle and half-hourly show the lowest values.Daily parametrisation reproduces the observations for average diurnal cycle better than the others in every occasion except the GPP for Sodankylä, where half-hourly tuning is better (verified by pointwise RMSE from the average diurnal cycle).We also notice that Sodankylä daily patterns, and to some extent Hyytiälä as well, are slightly out of phase.Our current understanding is that this is (at least partly) due to a slightly misaligned sensor (which can cause significant errors on high latitudes), measuring radiation fluxes.Fortunately this affects mainly the cost function for half-hourly tuning since it is the only one operating on the densest half-hourly timescale.

Dry event
Dry period in the summer 2006 can be clearly located by the massive drawdown in observed GPP, and to a lesser extent in ET, at Hyytiälä (Fig. 2).In a closer look at this event (Fig. 5) it is evident that none of our parametrisation schemes were able to capture it correctly.As it was with the time series, the overall structure of the daily time series during this event remains the same (there are no divergent aspects in the model output between the different tunings).
During the drought event (defined here as 31 July-15 August 2006), the soil moisture is on average 27 % lower for default, daily and half-hourly tuning and 40 % lower for seasonal tuning when compared to the corresponding values from other years -seasonal tuning has the lowest overall soil moisture.During this event the modelled soil moisture decreases monotonically for all tunings and reaches the lowest values on 13 August, after which it starts to rise.During the period the modelled ET and GPP are predominantly higher than the observations.WUE on the other hand follows the "observations" remarkably well and deviates from the observed only towards the end of the event when modelled ET drops to near-zero values, coinciding with the lowest modelled soil moisture values.Gao et al. (2016) examined deviation in the dependencies of GPP and ET to vapour pressure deficit (VPD) between model and observation results under the most severe soil moisture stress conditions at the end of the prolonged period of soil water scarcity (that occurred in 2006).This can be attributed to the lack of ex-  plicit dependence of the modelled stomatal conductance on the atmospheric humidity.

Conclusions
Initially we tuned the model to produce near-measured seasonal ET, GPP and especially maximum LAI to enhance the vegetation transpiration and to emphasise the response to precipitation.This was done successfully with seasonal tuning in the hopes of bringing forth the underlying model responses to dryness.With the consecutive daily and halfhourly tunings, we managed to improve the average diurnal cycles of both ET and GPP, but failed in reproducing the low ET and GPP levels during the dry event in 2006.Effectively we first (seasonal tuning) transferred water from soil moisture into (too high levels of) ET, and later (with daily and half-hourly tunings) transferred some of it back.
In addition to the maximum LAI ( max ) and maximum carboxylation rate (V C,max ), the most effective parameters in the seasonal tuning were the fraction of soil moisture above which transpiration is not affected by soil moisture stress (w tsp ) and the critical fraction of field capacity above which fast drainage occurs for soil water content (w dr ).The reduction in ET and GPP was mostly accounted for by lowering the approximate ratio of leaf internal CO 2 concentration to external CO 2 (f C3 ), which reduces the amount of carbon available for photosynthesis.For daily tuning ET was further reduced by the increase of the fraction of precipitation intercepted by canopy (p int ) and lower relative humidity fraction (w hum -air humidity is based on soil dryness).
Despite the fact that we were unable to enhance the dry response of the model, we are confident in saying that the algorithm itself worked well and performed as intended with the daily tuning providing the most reduction in model-data mismatch.We optimised 12 parameters simultaneously (with daily and half-hourly tunings) using eight fairly short chains (8000 samples).With daily tuning the resulting estimates are well matured, but with half-hourly tuning the parameter deviations are larger (which is probably due modelling inefficiencies and noise in measurements).Nevertheless, all optimisation procedures worked well with regard to what was optimised (seasonality, daily averages or diurnal cycle).
Recently, Knauer et al. (2015) found canopy conductance formulation to be a key factor in prescribing the transfer of carbon and water between terrestrial biosphere and the lower atmosphere.Additionally, Gao et al. (2016) found that during a prolonged period of soil water scarcity, the lack of explicit  dependence of the stomatal conductance on the atmospheric humidity is one of the contributing factors to this issue.Further studies into enhancing the dry response in JSBACH are needed and these studies should reflect these latest findings.

Data availability
The measurement data required to run and tune the model can be procured from the FLUXNET database (doi:10.17616/R36K9X).The JSBACH model is available to the scientific community under a version of the Max Planck Institute for Meteorology Software License Agreement (http: //www.mpimet.mpg.de/en/science/models/license/).For any questions regarding the simulations data, we encourage you to contact the author at jarmo.makela@fmi.fi.
drained (in addition to the limited drainage below this threshold).
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 = C h |v h | −1 ; these are heat transfer coefficient and horizontal velocity).
E ws = ρ q a − q s (T s , p s ) r a (A13) Transpiration from vegetation (T v ) is likewise formulated but additionally depends on the stomatal resistance of canopy (r).
T v = ρ q a − q s (T s , p s ) r a + r (A14) The stomatal resistance is given as a minimal stomatal resistance of the canopy without water stress (r min , depends on photosynthetically active radiation and LAI) divided by a water stress factor (f ws ).That is r = r min /f ws .The water stress factor depends on how much water is in the soil in relation to the maximum field capacity (w f = w s /w fc ) when compared to the limit when transpiration is no longer affected by soil moisture stress (w tsp ) and the permanent wilting point (w pwp ).
Evaporation from dry bare soil (E s ) is similarly defined as E s = ρ q a − hq s (T s , p s ) r a (A16) Here h is relative humidity at the surface relative to soil dryness: h = max w hum (1 − cos (π w f )) , min 1, q a q s (T s , p s ) . (A17) The total evapotranspiration is a weighted average of E ws , T v and E s , where the weights are based on, e.g., fill levels of reservoirs (similar to w f above) and vegetative fraction of the grid cell (veg max ).

Figure 2 .
Figure 2. Hyytiälä 7-day-running mean time series for different tunings for the first two summers of the validation period.Solid black line represents the observations.

Figure 3 .
Figure 3. Sodankylä 7-day-running mean time series for different tunings for the last two summers of the validation period.Solid black line represents the observations.

Figure 4 .
Figure 4. Average diurnal cycle from May to September for the validation period.

Figure 5 .
Figure 5. Daily averages for ET, GPP and WUE on a dry event in 2006 for Hyytiälä.

Table 1 .
Parameter descriptions with references to equations in Appendix A.
• C III LoGro phenology: memory loss parameter for calculating pseudo soil temperature.

Table 2 .
The highest correlations between parameters.

Table 4 .
Default and optimised parameter values using the last 20 000 samples (if no value is given, the parameter was not part of that tuning, and the default value was used instead).The percentage next to a parameter value is the effectiveness of that parameter for that tuning.The reference values for seasonal tuning are the default values and for daily and half-hourly tunings the seasonal values.