the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Optimal inverse estimation of ecosystem parameters from observations of carbon and energy fluxes
Debsunder Dutta
David S. Schimel
Ying Sun
Christiaan van der Tol
Christian Frankenberg
Canopy structural and leaf photosynthesis parameterizations such as maximum carboxylation capacity (V_{cmax}), slope of the Ball–Berry stomatal conductance model (BB_{slope}) and leaf area index (LAI) are crucial for modeling plant physiological processes and canopy radiative transfer. These parameters are large sources of uncertainty in predictions of carbon and water fluxes. In this study, we develop an optimal moving window nonlinear Bayesian inversion framework to use the Soil Canopy Observation Photochemistry and Energy fluxes (SCOPE) model for constraining V_{cmax}, BB_{slope} and LAI with observations of coupled carbon and energy fluxes and spectral reflectance from satellites. We adapted SCOPE to follow the biochemical implementation of the Community Land Model and applied the inversion framework for parameter retrievals of plant species that have both the C_{3} and C_{4} photosynthetic pathways across three ecosystems. We present comparative analysis of parameter retrievals using observations of (i) gross primary productivity (GPP) and latent energy (LE) fluxes and (ii) improvement in results when using flux observations along with reflectance. Our results demonstrate the applicability of the approach in terms of capturing the seasonal variability and posterior error reduction (40 %–90 %) of key ecosystem parameters. The optimized parameters capture the diurnal and seasonal variability in the GPP and LE fluxes well when compared to flux tower observations ($\mathrm{0.95}>{R}^{\mathrm{2}}>\mathrm{0.79}$). This study thus demonstrates the feasibility of parameter inversions using SCOPE, which can be easily adapted to incorporate additional data sources such as spectrally resolved reflectance and fluorescence and thermal emissions.
 Article
(11430 KB) 
Supplement
(8066 KB)  BibTeX
 EndNote
Terrestrial ecosystems play a very important role in regulating the carbon exchange over land surfaces (Falkowski et al., 2000; Schimel, 1995). Although they are known to be important sinks in buffering the increasing anthropogenic CO_{2} emissions (Friedlingstein et al., 2006; Sitch et al., 2015), there is a large variability and heterogeneity in the carbon exchange mechanisms, which are tightly correlated with interannual climatic variations (Cox et al., 2013; Liu et al., 2017). Moreover, terrestrial ecosystems also control the exchange of energy, water and momentum between the atmosphere and the land surface, thus regulating climate–ecosystem (carbon) feedbacks leading to amplification or dampening of regional and global climate change (Heimann and Reichstein, 2008). Measurements and modeling of carbon and water vapor fluxes over terrestrial ecosystems are therefore important to better understand these issues and account for the regional and global carbon and water budgets (Baldocchi et al., 1996, 2001; Sitch et al., 2003, 2008). Terrestrial ecosystem models have been used to study the carbon and water fluxes (Cramer et al., 2001; Kucharik et al., 2000; McGuire et al., 2001; Sitch et al., 2003); however, there are large uncertainties in fluxes associated with poorly quantified model parameters (Knorr and Heimann, 2001; Pappas et al., 2013; Rogers et al., 2017; Wramneby et al., 2008; Zaehle et al., 2005). Some of these parameters have temporal and spatial variability and are hard to measure directly over large scales (Dutta et al., 2017; Simioni et al., 2004; Wilson et al., 2000). For the majority of model implementations these parameters and their temperature dependence are represented as a single constant value according to plant functional types with little or no seasonal variability. In this study, we present an inversion approach which can be implemented with ecosystem models involving canopy physiological processes to better estimate the seasonal variability in photosynthesis and canopy structural parameters, which in turn can reduce the uncertainty in estimation of carbon and water fluxes over ecosystems.
The micrometeorological data from flux towers are extremely useful in understanding the biogeochemistry and thermodynamics of ecosystems (Baldocchi et al., 2000; McGuire et al., 2002). A number of approaches have been developed to model and estimate photosynthesis, respiration, energy balance, stomatal behavior, radiation transfer and turbulent gas exchange across the plant canopy on the basis of data from flux tower experiments (Oleson et al., 2010; Running and Coughlan, 1988; van der Tol et al., 2009). Detailed canopy models are often resolved into multiple layers, thus providing a better treatment of radiation regime and energy balance across the canopy (Dai et al., 2004; Wang and Leuning, 1998; van der Tol et al., 2009). At the heart of these models lies the leaflevel biochemical model of photosynthesis and carbon fixation (Farquhar et al., 1980) together with a stomatal conductance (most often the widely used Ball–Berry) model (Collatz et al., 1991b). The fluxes of carbon and water are tightly coupled through stomatal regulation and photosynthesis (Baldocchi, 1994; Collatz et al., 1992). Further, the processbased canopy models require some environmental drivers such as incoming shortwave and longwave radiation, air temperature, relative humidity, wind speed, and ambient CO_{2} concentration, along with a number of leaf and canopy parameters to simulate the fluxes of carbon in terms of gross primary production (GPP), flux of water or latent energy (LE), sensible heat (H), net radiation and others.
One of the most important ecosystem descriptors is the maximum rate of carboxylation (V_{cmax}), which is directly related to the concentration of the enzyme rubisco. V_{cmax} is a key parameter in the Michaelis–Menten kinetics for an enzymecatalyzed reaction of the substrates CO_{2} or O_{2} with ribulose1,5bisphosphate, representing the enzymelimited photosynthesis rate (Farquhar et al., 1980). Other ratelimiting photosynthesis parameters such as maximum electron transport rate (J_{max}) are generally parameterized with respect to V_{cmax}. The Ball–Berry equation calculates the stomatal conductance (g_{s}) for water vapor as a function of net assimilation, relative humidity, leaf surface CO_{2} concentration, minimum conductance and a proportionality constant called the Ball–Berry slope (BB_{slope}) (Ball et al., 1987; Beerling and Quick, 1995; Tanaka et al., 2002; Wullschleger, 1993). The BB_{slope} plays a crucial role in regulating the stomatal conductance and water use efficiency, and thus the surface energy fluxes in terms of partitioning the turbulent energy into LE and H fluxes. Thus, it is a crucial parameter regulating the tradeoff between carbon gain and water loss, e.g. during drought conditions (Monteith and Unsworth, 2007). The leaf area index (LAI) is a canopy structural and key ecosystem variable in most terrestrial biosphere models, which determines interception of radiation as well as photosynthesis and energy exchange across the canopy (Chen et al., 1997). The parameters V_{cmax} and BB_{slope} can be determined experimentally from leaflevel gas exchange measurements and generated A−C_{i} curves (Tanaka et al., 2002; Wullschleger, 1993; Xu and Baldocchi, 2003). LAI can be estimated from destructive and nondestructive optical methods (Chen et al., 1997; Dutta et al., 2017; Myneni et al., 1997), as well as inversion approaches on spectrally resolved reflectance data from satellite and airborne platforms (Houborg et al., 2007; Jacquemoud et al., 1995). However, these measurements are much more complex and labor intensive, being measured less frequently than flux tower observations.
Inversion of detailed processbased models using observations of carbon and energy fluxes could thus yield these key ecosystem parameters. Processbased models such as the Soil Canopy Observation, Photochemistry and Energy fluxes (SCOPE) model (van der Tol et al., 2009) can simulate the radiative transfer and the fluxes of carbon and energy vertically resolved within the canopy. Our hypothesis is that the inversion of a detailed vertically resolved canopy model such as SCOPE with multiple layers consisting of sunlit and shaded fractions together with fully spectrally resolved radiation regime and energy balance computations (van der Tol et al., 2009) is able to retrieve the ecosystem parameters accurately using observations of carbon and energy fluxes, and in the future remote sensing data, as SCOPE can model the spectrally resolved shortwave reflectance, thermal emission and solarinduced chlorophyll fluorescence.
A few studies have used inversion approaches to extract ecosystem parameters from flux (Reichstein et al., 2003; Schulze et al., 1994) and reflectance (Quaife et al., 2008) measurements but not yet to constrain all three key parameters (V_{cmax}, BB_{slope} and LAI) simultaneously using the fluxes of water and carbon. A previous study by Wolf et al. (2006) used a deterministic linear leastsquares inversion method to estimate the key ecosystem parameters (V_{cmax}, BB_{slope}, LAI and respiration rate) using the net ecosystem exchange (NEE) and sensible and latent heat fluxes. The approach assumed a simple model of radiationdriven photosynthesis, respiration and energy balance using a twocomponent (sunlit and shaded) canopy. The optimization used total energy (H + LE) to fit LAI values, the NEE to fit V_{cmax} and respiration rate and energy difference (H − LE) to fit BB_{slope}. In comparison to a deterministic approach, stochastic Monte Carlo approaches (Knorr and Kattge, 2005; Mackay et al., 2012; Ricciuto et al., 2008; Xu et al., 2006) constrain a number of parameters (including the photosynthetic parameters) using eddy covariance observations but assuming them to be time invariant. These studies consider multiple temporal resolutions such as seasonal or halfyearly resolutions and present a range of parameters without providing a definite error characterization as Bayesian methods (Wu et al., 2009). Moreover, since the stochastic methods sample the probability distribution in parameter space, they are better suited to nonlinear models, but the associated computational costs can often be prohibitive.
In this study, we develop an inversion framework for estimating the temporal dynamics of key ecosystem parameters using the SCOPE model, which represents detailed plant physiological processes including suninduced chlorophyll fluorescence (SIF). SIF is chlorophyll reemission during photosynthesis, acts as a direct probe into photosynthesis measurable from space and is strongly correlated with fluxbased GPP estimates at canopy to ecosystem scales (Flexas et al., 2002; Frankenberg et al., 2011). Thus, the SCOPEbased inversion approach has the flexibility and advantage of incorporating towerbased observations of fluxes including SIF as well as spectrally resolved reflectance and thermal emissions for optimal estimation of a wide range of ecosystem parameters. In this paper, we first focus on the conceptual framework of parameter inversion using SCOPE followed by parameter retrieval examples, with specific objectives as follows:

implementation of photosynthesis model and its temperature dependencies consistent with a wellaccepted major Earth system model (Community Land Model CLM4.5) in SCOPE;

development of a Bayesian nonlinear inversion framework using SCOPE to estimate ecosystem parameters using eddy covariance flux observations;

demonstrating the retrieval and posterior error reduction of key ecosystem parameters using (i) observations of carbon and water fluxes and (ii) combining flux observations together with satellite reflectance across different ecosystems.
The rest of the paper is organized as follows. Section 2 provides a brief overview of the SCOPE model and the new implementation of photosynthesis and its temperature dependencies. Section 2.2 provides a comparison of the old and new photosynthesis implementations in SCOPE. Sections 3, 4 and 5 describes the formulation of the inverse problem followed by linearization of the forward model and mechanisms of the retrieval algorithm. Section 6 describes the results of the inversion framework across three different ecosystems and finally Sect. 7 provides a discussion summary and conclusions.
The SCOPE model (van der Tol et al., 2009) is an integrated 1D vertical radiative transfer and energy balance model. The model utilizes the spectrally resolved visible to thermal (0.4 to 50 µm) infrared irradiation at the canopy top to derive the fluxes of water, energy, carbon dioxide and vertical profiles of temperature as a function of canopy structure and weather variables. The four most important SCOPE modules represent (i) radiative transfer of incident solar radiation and generated fluorescence within the leaf (Fluspect), (ii) radiative transfer of incident direct and indirect solar radiation (0.4–50 µm), (iii) radiative transfer of internally generated thermal radiation by vegetation and soil (Verhoef et al., 2007), (iv) an energy balance module and (v) a radiative transfer module for computing the topofcanopy radiance spectrum of fluorescence from leaflevel chlorophyll fluorescence. SCOPE resolves topofcanopy incoming and outgoing shortwave radiation and reflectance in the spectral range of 400 to 2500 nm at 1 nm wavelength bands. Further, it also computes the spectrally resolved fluorescence emission in the range of 650 to 850 nm at 1 nm wavelength bands.
One important aspect is that SCOPE relaxes the assumption of constant temperatures for the sunlit and shaded fractions of the leaves across the different canopy layers. This is true when we consider different orientations and their vertical positions in the canopy. Therefore, an iterative solution scheme is implemented in SCOPE as stomatal conductance affects leaf temperature, which in turn affects photosynthesis (and thus again stomatal conductance). Thus, the fully integrated thermal radiative transfer and energy balance modules allow feedback between leaf temperatures, photosynthesis, chlorophyll fluorescence and radiative fluxes.
2.1 The SCOPE biochemical module
The SCOPE biochemical module is a submodule of the energy balance routine, which provides an iterative solution of the photosynthesis, energy balance, net radiation and heterogeneous skin temperatures for a particular net external forcing. The main functions of the biochemical module include leaf temperature dependent computation of photosynthesis and fluorescence. Some of the photosynthesis parameterizations in the current version of the SCOPE model (V1.70) are outdated and more in line with previous versions of the Community Land Model (CLM version 4) or based on a mix of other model implementations. CLM is a communitydeveloped land model which focuses on the modeling of land surface processes including biogeophysics, carbon cycle, vegetation dynamics and river routing. Specifically, the main modifications in the more recent CLM (version 4.5, CLM4.5) (Lawrence et al., 2011; Oleson et al., 2013) include updates to the canopy radiation scheme and canopy scaling of leaf processes, colimitations on photosynthesis, revisions to photosynthetic parameters (Bonan et al., 2011), temperature acclimation of photosynthesis, and improved stability of the iterative solution in the photosynthesis and stomatal conductance model (Sun et al., 2012). CLM4.5 implements a multilayer canopy modeling framework with coupled photosynthesis (Farquhar et al., 1980) and Ball–Berry stomatal conductance models similar to the SCOPE framework.
The main inconsistencies of SCOPE (V1.70) with the CLM4.5 parameterizations are as follows:

Similar, generic temperature response functions are implemented for both C_{3} and C_{4} species, excepting V_{cmax}, and further it uses a Q_{10}based exponential function with the same functional parameters for computing the temperature response of the various photosynthetic parameters.

There is no J_{max} (maximum potential electron transport rate, ETR) or its temperature dependence in the computation of light limited C_{3} photosynthesis rate.

The net assimilation, internal CO_{2} concentration and stomatal conductance (AC_{i}g_{s}) iterative solution method is not quite robust or was lacking in the previous versions, with the V1.70 implementation being complicated and unpublished.
We therefore attempt to improve the SCOPE biochemical module by implementing the photosynthesis and temperature dependence of the photosynthetic parameters according to the well established and widely used CLM4.5. All the detailed implementation steps and equations for modeling the photosynthesis and temperature dependence primarily as per Bonan et al. (2011) is presented in detail in the Appendix A and B. Within the inverse framework described later, we only invert V_{cmax} at the reference temperature of 25^{∘} and apply the given temperature dependencies. Any systematic difference in the temperature function could thus alias into the derived V_{cmax}. Overall, the major new updates made to the model (biochemical module) are as follows:

The electron limited photosynthesis rate A_{j} is computed using the potential ETR J, which is obtained by solving the smaller root of Eq. (A5), comprising the light utilized in photosystem II (I_{PSII}) and the maximum potential ETR (J_{max}).

The light limited photosynthesis rate for C_{4} is given by Eq. (A3).

The temperature dependence of photosynthetic parameters (Bonan et al., 2011) now uses the activation, deactivation energies and entropy terms in the temperature response and high temperature inhibition functions (Leuning, 2002) (see Appendix B for details). The temperature response of C_{3} (Bernacchi et al., 2001; Leuning, 2002) and C_{4} photosynthesis is represented by Eqs. (B1)–(B5).

Finally we also incorporate a new simplified implementation of AC_{i}g_{s} iterations (Sun et al., 2012) and include the computation of oxidative photosynthesis (Bernacchi et al., 2001) within the photosynthesis model. See Appendix B1 for details.
In the following section, we demonstrate the photosynthesis results with the newer photosynthesis and temperature dependence implementation as well as its comparison with the previous version (V1.70) in SCOPE for different ecosystems.
2.2 Comparison of current and previous photosynthesis implementations in SCOPE
Figure 1 shows the temperature response functions for V_{cmax} for both C_{3} (left) and C_{4} (right) photosynthetic pathways. The functions of mean temperature response, high temperature inhibition and the 1σ variance as per the different photosynthesispathwaydependent parameterization (e.g. activation energy, deactivation energy, entropy) are shown according to Leuning (2002). The new temperature dependency parameterizations follow the temperature functions and high temperature inhibition for C_{3} and the Q_{10} functions for the C_{4} pathways. We have also shown the temperature dependence of V_{cmax} from the previous (old) implementation of the SCOPE model (V1.70). The differences in the net response at both lower and higher than optimal temperature can be clearly identified in the figure for both C_{3} and C_{4} species. It can be observed that the difference in temperature response is more for C_{4}, and clearly the maximum is in the leaf temperature range 30–40^{∘}C; however, it continues into the higher temperatures as well. Moreover, it can be noted that the overall shapes of the response functions are nearly identical (with some lag) for the different parameters for the previous SCOPE implementation compared to the newer implementation as per CLM4.5 (Bonan et al., 2011).
A number of analyses were performed to study the differences in net response of canopylevel carbon and energy fluxes for both C_{3} and C_{4} species from the SCOPE model due to modification in photosynthesis implementation (old and new) and its temperature dependence (see Supplement Sect. S1).
Figure 2 shows the comparison between the old and new SCOPE versions as the ratio of overall canopy GPP, which is defined as ${f}_{\mathrm{GPP}}=\frac{{\mathrm{GPP}}_{\mathrm{new}}}{{\mathrm{GPP}}_{\mathrm{old}}}$. This ratio is further represented as a function of the three most important forcing variables: PAR (photosynthetically active radiation), canopy temperature and VPD (vapor pressure deficit). The results for the Missouri Ozark site (see Sect. 6.4.1 for site details) with C_{3} plant species for the year 2009 are presented in Fig. 2. For this analysis, the SCOPE model simulations are computed for the entire growing season and the f_{GPP} values are binned according to the PAR temperature (for specific VPD ranges) and PAR VPD (for specific temperature ranges) as 2D histograms, of which only the mean (${f}_{{\mathrm{GPP}}_{\mathrm{mean}}}$) is represented in Fig. 2.
We find that over the larger parts of the domain of random variables, f_{GPP} is around 1 and the maximum change in overall GPP is around 25 %. From Fig. 2, it can be observed that in the case of C_{3} species, for the combinations of higher canopy temperature and low VPD values (panel a), the new GPP values remain the same or are reduced by about 5 %, although from Fig. 1 we find an increase in the V_{cmax} response at >25 ^{∘}C temperature, which may indicate photosynthesis being limited by light instead of the enzyme rubisco. For the combination of low canopy temperature and lower VPD values (panel c), f_{GPP} values are close to 1 (except for very low PAR/VPD values and with a maximum of about 4 %–6 % increase), which can be explained by an almost identical V_{cmax} response at lower temperatures in Fig. 1. At high PAR values with higher temperatures (25–30 ^{∘}C) and low VPD values (panel d) we find that GPP increases by about 6 %–14 %, which can be directly explained by the new increased V_{cmax} response in that temperature range as indicated in Fig. 1. The results for a similar comparative study with C_{4} species is presented in the Sect. S1. Overall, we find that the new model implementation of photosynthesis and its temperature dependence as well as A−C_{i} iterations works well and only result in moderate, yet noticeable changes. It also underlines that tabulated model parameters can only be optimized for a specific model implementation, which is not necessarily universally transferable to other carbon cycle models.
The problem of ecosystem flux computation (e.g. GPP, latent energy) from meteorological variables (e.g. VPD, air temperature, relative humidity) and other ecosystem parameters can be represented as follows:
where $\mathcal{F}\left(\right):{\mathit{X}}^{\mathbf{\prime}}\u27f6\mathit{Y}$ is a functional representation of the model, which maps the model input and parameter space (X^{′}) quantitatively to the space of ecosystem fluxes (Y), and ε represents the residual error which includes the precision error, the model error and random errors. In our case, SCOPE represents the forward model ℱ(), which is complex and moderately nonlinear, representing a range of physics and canopy physiological processes. We can further represent our forward problem as follows:
where X represents the state vector of parameters to be retrieved, p ($\mathit{X},\mathit{p}\subset {\mathit{X}}^{\mathbf{\prime}}$ and ${\mathit{X}}^{\mathbf{\prime}}=\mathit{X}\cup \mathit{p}$) is a vector of parameters which represents those quantities that influence the measurement and are known with some accuracy but are not meant to be retrieved. We call these parameters the forward functional parameters. In our example p represents the set of all fixed model (SCOPE) parameters not involved in the retrieval. The error term ε represents the measurement noise (e.g. noise or errors in the flux measurements). Given a set of measurements Y, the optimal state vector $\widehat{\mathit{X}}$ can be obtained by a generalized inverse method ℛ represented as follows:
where $\widehat{\mathit{p}}$ represents the best estimate of the forward function parameters. The parameters X_{a} and c represent the parameters that do not appear in the forward function but they do affect the retrieval and are associated with uncertainties. X_{a} represents the prior estimate of X and c represents any other parameters in the retrieval scheme as a catchall for anything else that is used in the retrieval method, which also includes the convergence criteria.
A basic prerequisite for inverting the forward model is to compute its sensitivity with respect to input parameters, i.e. the partial derivatives with respect to all the state vector elements (Jacobi matrix). For linear models, the Jacobians are independent of the actual state. In our case, the SCOPE forward model is moderately nonlinear and its Jacobians need to be computed numerically as analytical methods are currently lacking and hard to implement given some peculiarities in the FvCB equations.
With the Jacobian matrix and a simple forward model call, we can thus write a firstorder Taylor expansion for the forward model
where X_{l} is an arbitrary linearization point, and $\frac{\mathit{\delta}F}{\mathit{\delta}\mathit{X}}$ is the partial derivative or Jacobian at the point X=X_{l}.
In the remainder of the paper, we will omit the vector of forward model parameters p which are not a part of the retrieval framework. For the nonlinear problem we use the maximum a posteriori approach. The Bayesian solution for the nonlinear inverse problem where the forward model is a general function of the state, the measurement error is Gaussian (S_{ε}) and with a prior estimate of the state (X_{a}) with a Gaussian uncertainty in the prior state (S_{a}) (Rodgers, 2000) can be represented as follows:
where c^{′} is a constant. Our aim is to find the best estimate of the state vector $\widehat{\mathit{X}}$ (denoted as X henceforth) and an error characterization that describes the posterior PDF. The Gauss–Newton iteration steps for determining the state vector is given by the following:
where K_{i} is the Jacobi matrix. A brief derivation of Eq. (6) is presented in Appendix C; for a more indepth treatment the reader is referred to Rodgers (2000).
5.1 Levenberg–Marquardt method
In general, the Gauss–Newton iterations discussed previously finds the minimum in one step if the cost function is quadratic with respect to X. However, in our case the cost function is not perfectly quadratic and the initial guess potentially far away from the solution, thus requiring multiple iterations. In addition, the nonlinearity of the problem sometimes results in steps that would actually increase rather than decrease the fit quality. In order to overcome this issue Levenberg (1944) and Marquardt (1963) proposed the following iteration for nonlinear leastsquares problem:
where γ_{i} is chosen at each step to minimize the cost function and D is a diagonal scaling matrix to scale the elements of the state vector. It can be noted that for γ_{i}→0, Eq. (7) leads to a Gauss–Newton iteration step and for γ_{i}→∞ Eq. (7) tends to steepest descent and further the step size tends to 0. It is also expected that the cost function will decrease corresponding to the decrease in γ_{i} from infinity to zero. The value of γ_{i} is sequentially updated at each iteration by evaluating the change in cost function. Here, we follow the general recommendations as outlined in Marquardt (1963) and Rodgers (2000).
The guidance for choosing the scaling matrix D is that it must be positive definite. For the current problem we choose it to be ${\mathbf{S}}_{\mathrm{a}}^{\mathrm{1}}$ (as in Rodgers, 2000) and apply the Levenberg–Marquardt (LM) modification to the Gauss–Newton method (iteration Eq. C8), resulting in the following iterative inversion scheme:
5.2 A moving window setup of the inversion problem using flux tower observations
The top row of Fig. 3 shows the SCOPE model simulations of GPP, LE, H and SIF for one day (3 August 2010) in the growing season for C_{4} corn using data from the Nebraska Mead1 flux tower site with parameter values V_{cmax}=50 µmols m^{−2} s^{−1}, BB_{slope}=7 and LAI =4. The second, third and fourth rows from the top show the numerically computed partial derivatives of GPP, LE, H and SIF with respect to the parameters using SCOPE with positive perturbations δV_{cmax}=5 µmols m^{−2} s^{−1}, δBB_{slope}=1 and δLAI =0.5. Each column of Fig. 3 represents a row of Jacobian matrices used for the inversions. The figure clearly demonstrates the influence of each of the parameter variables in the state vector (X) on the modeled fluxes (F(X)). We can observe the counteracting nature of variables and the fluxes from the Jacobian. For example, for LE flux, BB_{slope} has a positive gradient but LAI has a negative gradient. The decrease in LE is attributed to less radiation reaching the soil and a corresponding increase in soil aerodynamic resistance. In this case the canopy resistance goes up but does not compensate for the decreased soil evaporation and results in low sensitivities. Similarly we find V_{cmax} has a positive gradient for GPP but a negative one for LE, which may again be attributed to the soil evaporation responding to soil temperature. It can be noted that the nature of these sensitivities at the canopy level are sometimes counterintuitive from their leaflevel mechanisms and may vary depending on environmental conditions, such as incoming PAR as well as air temperature and vapor pressure deficit. This also creates diversity in the Jacobians over the diurnal cycle, which allows us to derive more than two parameters from two sets of measurements (GPP and LE). In Fig. 3, we have not only shown derivatives of GPP and LE but also H and SIF (not used here). In this paper, we outline the general framework of parameter inversion, which can easily be modified to make use of more measurements such as H, SIF, reflectance or thermal emissions, all of which can be modeled with SCOPE.
For setting up the observation vector Y_{m×1} (see Eq. 8), we use observations of carbon and latent energy fluxes from eddy covariance tower time series records. The observational error matrix (S_{ε[m×m]}) is assumed to be a diagonal matrix and computed using noise standard deviation as 10 % of the halfhourly to hourly observations (Leuning et al., 2012). We use an initial prior state value of the state vector (X_{a[n×1]}) as well as the prior error covariance matrix (S_{a[n×n]}). As mentioned previously, the Jacobian matrix K_{m×n} is computed numerically by a small perturbation to the value of the state vector X_{i}+δ (see Fig. 3) at a particular iteration step. The observed (Y) and modeled (F(X)) fluxes in the inversion framework are set up as a long concatenated vector as shown in Fig. 4. The concatenation of different flux variables is done using a time filter to represent the part of the day we wish to include in the retrieval framework as illustrated in Fig. 4. This is logical as we have already demonstrated in Fig. 3 that the gradients are variable throughout the day. Ideally, the time filter applied for concatenating the data should capture the maxima and a range of variations in the gradients, but at the same time reduce the data points to make the retrieval computationally efficient and further tend towards providing stable solutions (retrievals) of the parameter values. Further, the time filter helps to eliminate the nighttime anomalies in the observations for accurate parameter estimation. For other observations such as spectral reflectance a daily noontime average is suitable for concatenation in the observations Y. The assumptions behind the longterm (seasonal) retrieval of important ecosystem and plant physiological parameters is that these parameters change significantly over the growing season but at a slower rate compared to and in response to the environmental and meteorological forcing. Thus, the ecosystem parameters can be assumed to be constant over some finite time window. We implement this assumption to set up our inverse parameter retrieval framework for finite Nday contiguous moving windows over the entire growing season (Fig. 4). We extend the 1day diurnal setup of Y, F(X) and K as shown in Fig. 3 to multiple days for setting up the Nday windows as illustrated by color coding in Fig. 4. After computing the necessary vectors and matrices for the Nday window, iterations are performed by applying the LM algorithm until convergence to obtain the posterior estimation of the state vector. The retrieval window is moved over to the contiguous next N days and the process is repeated. The retrieval proceeds thus for the entire length of the growing season (Fig. 4). For our retrieval example, we choose a 3day moving window which seems optimal for the plant response in terms of the photosynthesis parameters (V_{cmax}, BB_{slope} and LAI) towards the change in environmental drivers.
5.3 Error characterization and convergence criteria for the retrievals
As mentioned in Sect. 5.1, we have selected a convergence criteria for the parameter retrievals in each of the moving windows based on the ratio of the true error to the expected error for each of the iteration steps. The total error minimized for the retrieval is given by $[\mathit{Y}F({\mathit{X}}_{i}\left){]}^{T}{\mathbf{S}}_{\mathit{\epsilon}}^{\mathrm{1}}\right[\mathit{Y}F\left({\mathit{X}}_{i}\right)]+[\mathit{X}{\mathit{X}}_{\mathrm{a}}{]}^{T}{\mathbf{S}}_{\mathrm{a}}^{\mathrm{1}}[\mathit{X}{\mathit{X}}_{\mathrm{a}}]$. However, for testing the convergence within each iteration step, we use the method suggested by Rodgers (2000), which adapts the Levenberg–Marquardt parameter depending on the nonlinearity of the forward model.
After convergence, the posterior error covariance matrix for the retrieved state vector $\widehat{\mathit{X}}$ can be computed as follows:
The reduction in error is defined as follows:
where S_{jj} and ${\mathbf{S}}_{{\mathrm{a}}_{jj}}$ represent the diagonal elements in the posterior and prior error covariance matrices respectively. In our LM retrieval process, we use the retrieved state vector $\widehat{\mathit{X}}$ of the previous window as the first guess (but not prior) for the current window. This saves computational cost and is based on the assumption that our state vector varies smoothly in time.
In this section, we discuss the results of optimal parameter estimation by applying the Bayesian inversion framework to three different ecosystems. The aim is to demonstrate the applicability for the retrieval (as well as capturing the seasonal variability) of canopy structural and photosynthesis parameters using carbon and water fluxes, and to further compare and contrast the results across the different ecosystems. In order to demonstrate the greater potential of SCOPE in modeling spectrally resolved reflectance (not found in other general carbon cycle models) and versatility of the inversion framework we have also incorporated Moderate Resolution Imaging Spectrometer (MODIS) satellite reflectance bands in the retrievals. We further demonstrate how reflectance and fluxes are able to better constrain parameters such as V_{cmax}, BB_{slope} and LAI compared to just using flux tower observations.
6.1 Data filtering criteria in the moving window retrievals
Apart from the overall algorithmic steps as described previously, we apply the following filter criteria on the results and the data for a computationally efficient retrieval.

In constructing the observation vector Y we apply a timeoftheday filter (e.g. data between 09:00 and 16:00 LT and so on) for the initial forward SCOPE model.

For computing the Jacobians, a PARbased threshold (PAR >100 µmols m^{−2} s^{−1}) is applied to ensure sensitivity of the measurement vector with respect to state vector variations and to minimize the occurrence of unreasonable flux tower data (high fractional errors).

A filter is implemented to check and ensure that the state vector remains positive at every iteration. If somehow due to a small enough γ the state vector is negative, the γ value is adjusted in an iterative manner to keep it within bounds.
6.2 MODIS satellite reflectance data
We use the daily MODIS MCD43A reflectance product in this study (Schaaf and Wang, 2015). The spatial resolution of the dataset is 500 m and bands 1 and 2 (red and NIR) centered at 620 and 841 nm respectively were used in the inversion. These data are adjusted using a bidirectional reflectance distribution function to model the values as if they were collected from a nadir view. Figure 5 shows the distinct seasonality in greenness which is represented by NDVI over the two sites. SCOPE models the full Nadir VSWIR spectral reflectance (400 to 2500 nm) from which values corresponding to the two MODIS reflectance bands are extracted and used concurrently with the observations in the inversion framework.
6.3 Retrieval results for the Nebraska Mead1 site
6.3.1 Site description
The Nebraska Mead1 site is a part of the AmeriFlux network, located in Lincoln, Nebraska, and is one of the three cropland sites at the University of Nebraska Agricultural Research and Development Center, with continuous data records from 2001 until the present (Suyker et al., 2005). This site is a continuously irrigated corn (C_{4}, species) crop site, with mean annual precipitation of 790 mm and mean annual temperature of 10.07 ^{∘}C. We choose the year 2010 and an hourly time resolution for the analysis. The site meteorology and forcing variables relevant to the SCOPE inversion retrievals are shown in Fig. 6. The top two panels show the environmental forcing variables which are used as input (except precipitation) in the SCOPE simulations. The bottom panel represents carbon (GPP) and energy (LE, H) fluxes, which are used to construct the observation vector Y. The figure indicates that the growing season extends from around June through September, coinciding with high temperature, VPD and net radiation. We focus on the retrieval of the parameters V_{cmax}, BB_{slope} and LAI during this entire growing season.
6.3.2 Inversion parameters and results
For each of the retrieval windows, the prior value of the state vector along with prior errors and daytime duration, which is used for filtering the GPP and LE observations, are shown in Table 1. Here, we use a purely diagonal prior error covariance matrix, with zero offdiagonal elements. Figure 7 shows the retrievals of parameters V_{cmax}, BB_{slope} and LAI. The grey time series of GPP and LE values in the background are the actual filtered values used for constructing the observation (Y) vector corresponding to each retrieval window. The dashed lines indicate the retrieved parameters using only GPP and LE fluxes. The solid lines indicate the retrieved parameters using the fluxes together with MODIS reflectance. The orange (with fluxes only) and brown (with fluxes and MODIS red and NIR reflectance) lines show the result of posterior simulations of fluxes with the optimized parameters. These lines represent the absolute daily average posterior simulation errors ($\mathrm{\Delta}=\mathrm{Observed}\mathrm{Posterior}$) in the fluxes without and with the use of MODIS reflectance along with the flux observations in the inversions.
We find a seasonal variability in the retrieved parameters, which follow a similar pattern in GPP or LE. In particular, the retrieved LAI as well as V_{cmax} shows a similar seasonality to GPP. There is also some variability in retrieved BB_{slope}, which is correlated with LE observations. We found that including MODIS reflectance places better constraints on the parameters during the peak of the growing season, with much less variability in retrieved LAI and V_{cmax}.
As expected, the optimized parameters using just flux observations (dashed lines) are quite sensitive to the variation in GPP and LE, for example around DOY 190 and 210, where there is sudden dip in the retrieved V_{cmax}. Including the MODIS reflectance (solid lines) in the inversions alleviates most of these large variabilities due to fluctuations in the observed fluxes or meteorology. Moreover, with this additional MODIS reflectance constraint the range of variability for all the three parameters V_{cmax}, BB_{slope} and LAI is more realistic. The comparison of posterior simulations (brown and orange lines) indicates the net errors in the prediction of fluxes are quite similar in both cases (with and without MODIS reflectance) with the difference between the two being δ_{Δ}≤15 % during the middle of the growing season. This may indicate that in this example there is an equifinality in the posterior simulation of fluxes with the retrieved parameters, which gets alleviated with the reflectance data. We find that the MODIS reflectance better constrains LAI during the beginning of the growing season between DOY 160 and 180. The unexpectedly large increases in BB_{slope} and LAI around DOY 250–260 may be partially attributed to the largest rainfall events (see Fig. 6). Part of this variability and correlation between BB_{slope} and LAI may also be due to the diminishing role of soil evaporation (parameterized by a single resistance in SCOPE) with increasing LAI. Another part may be due to evaporation from the wet canopies, which is not currently represented in SCOPE. This may cause the inversion to overestimate BB_{slope}, even though it would not represent the gas exchange through the stomata. From the inversion results it is clear that all three parameters V_{cmax}, BB_{slope} and LAI are much better constrained (with more realistic values and better seasonal variability) with the assimilation of reflectance data together with fluxes in the optimal estimation framework.
Figure 8a shows the final posterior error reduction (ζ_{i}Eq. 10) of the retrieval iterations for each moving window. The dashed lines indicate retrieval results using GPP and LE observations only and the solid lines show retrievals that use reflectance observations in addition. The value of ζ_{i} is computed from the diagonal elements of the posterior error covariance matrix. We find a significant reduction in the posterior errors of the variables in the state vector. There is a strong seasonality in ${\mathit{\zeta}}_{{V}_{\mathrm{cmax}}}$ and ζ_{LAI} values and moderate to no seasonality for the ${\mathit{\zeta}}_{{\mathrm{BB}}_{\mathrm{slope}}}$. It can be clearly seen that the posterior error reduction is significantly greater (∼50 %) when combining the reflectance data with the flux observations. The error reduction provides more confidence in the retrieved parameters, which are also more realistic. The posterior error covariance matrix also indicates whether the retrieved parameters are truly independent (as in the case of a diagonal matrix) or whether they covary (indicated by significant offdiagonal elements). The error correlation is given by ${\mathit{\rho}}_{x,y}=\frac{\mathrm{COV}(x,y)}{{\mathit{\sigma}}_{x}{\mathit{\sigma}}_{y}}$ and should be considered when interpreting covariations of retrieved parameters as the nature and magnitude of the associations between the variable pairs are true for the retrievals only and may or may not represent the behavior of the variable pairs in nature due to different environmental conditions. Figure 8b shows the growing season error correlation patterns between the three parameters from the retrievals. From the results which include reflectance data with flux observations, it can be seen that during the peak growing season ${\mathit{\rho}}_{{\mathrm{BB}}_{\mathrm{slope},\mathrm{LAI}}}$ and ${\mathit{\rho}}_{{V}_{\mathrm{cmax},\mathrm{LAI}}}$ have opposite signs, indicating slightly positive and negative association between the variable pairs in the retrievals. These trends are also true when only flux data are used for the retrievals. However, in comparison ${\mathit{\rho}}_{{V}_{\mathrm{cmax}},{\mathrm{BB}}_{\mathrm{slope}}}$ is mostly zero but has opposite signs with and without reflectance data and thus indicates both favorable and competing effects during the middle of the growing season.
Finally, Fig. 9 demonstrates the net improvement in canopy GPP and LE fluxes due to the optimized state vector using flux observations only over their prior values. The first column represents the diurnal and seasonal variability in the time series of GPP and LE fluxes with optimized and unoptimized parameters and further its comparisons with flux tower values. The right column represents the onetoone comparisons of the same. We find a significant improvement in the estimation of GPP (R^{2}=0.94) and the optimized parameters are able to capture the growing seasonal variability well as measured by flux tower observations (slope =1.04). The improvement in modeling the fluxes with the posterior over the prior value of the state vector is also captured by the χ^{2} error statistic (${\mathit{\chi}}^{\mathrm{2}}={\sum}_{i=\mathrm{0}}^{k}\left(\right({y}_{\mathrm{i}}F\left({x}_{\mathrm{i}}\right))/{\mathit{\sigma}}_{\mathrm{i}}{)}^{\mathrm{2}}$) for the prior (unoptimized) and posterior (optimized) simulations. The corresponding values are ${\mathit{\chi}}_{\mathrm{GPP}\mathrm{opt}}^{\mathrm{2}}=\mathrm{9382}$, ${\mathit{\chi}}_{\mathrm{GPP}\mathrm{unopt}}^{\mathrm{2}}=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{728}$, ${\mathit{\chi}}_{\mathrm{LE}\mathrm{opt}}^{\mathrm{2}}=\mathrm{19}\phantom{\rule{0.125em}{0ex}}\mathrm{235}$ and ${\mathit{\chi}}_{\mathrm{LE}\mathrm{unopt}}^{\mathrm{2}}=\mathrm{20}\phantom{\rule{0.125em}{0ex}}\mathrm{554}$.
6.4 Retrieval results for the Missouri Ozark site
6.4.1 Site description
The Missouri Ozark site is also a part of the AmeriFlux network and is located in the University of Missouri Baskett Wildlife Research area, situated in the Ozark region of central Missouri. It is uniquely located in the ecologically important transitional zone between the central hardwood region and the central grassland region of the US (Gu et al., 2006). This site has a mean annual precipitation of 986 mm and a mean annual temperature of 12.11 ^{∘}C and has continuous data records from 2004 until the present. It is a deciduous broadleaf forest site comprised of C_{3} plant species. We use halfhourly datasets from the year 2007 and 2009 in the present analysis. The site meteorology and forcing variables relevant to the SCOPE inversion retrievals are shown in Figs. S4 and S5 in the Supplement. From the meteorological data it can be seen that the year 2009 was a normal wet year and the year 2007 was a year with a midsummer drought around DOY 250. This is also reflected in observed GPP and LE fluxes with two distinct peaks in the growing season, caused by a latesummer drought and associated low productivity around DOY 250. This decrease in productivity is not distinguishable from the MODIS reflectance data in Fig. 5, which indicates that plants maintain greenness during this time, making this a unique test case for our inversion setup as the V_{cmax} fits reflect the stress factor as well, which is usually applied to downscale the physiological V_{cmax} during environmental stress. We focus on the retrieval of the parameters V_{cmax}, BB_{slope} and LAI during this entire (longer) growing season. For both years we demonstrate the parameter retrievals using GPP and LE fluxes only as well as compare our retrievals using additional constraints of MODIS red and NIR reflectances (Fig. 5).
6.4.2 Inversion parameters and results for the year 2009
The assumed prior values of the state vector, prior errors and daytime duration which are used for filtering the GPP and LE observations for the retrieval windows are shown in Table 1. Compared to Mead, the retrieval for this site is carried out over a much longer duration, covering almost the entire year. Figure 10 (like Mead1, Fig. 7) shows the comparison of parameter retrievals V_{cmax}, BB_{slope} and LAI using (i) only flux observations and (ii) flux observations combined with two MODIS reflectance bands. We find an overall strong seasonality with realistic values of V_{cmax}, BB_{slope} and LAI, following the patterns in GPP and LE. The beginning (and end) of the growing season shows increasing (decreasing) trends in V_{cmax} and LAI retrievals, which coincide well with the increase (decrease) in MODIS NDVI observations (Fig. 5). The retrieved LAI variability helps to explain the rapid appearance of new leaves in the spring (March–April) and their disappearance around fall (October–November). The MODIS reflectance data help to better constrain LAI and V_{cmax}, which is specifically evident around DOYs 140–150 and 200–250. The sharp increase in retrieved V_{cmax} when using just the flux observations around DOY 148 may be attributed to the sharp peaks in GPP; however, the reflectance observations clearly help to better constrain the state vector. There are some issues with the retrieval for windows around DOY 175 and 275 (BB_{slope}); this may again correspond to the large precipitation (see Fig. 4) events (e.g. overestimation of BB_{slope}) around these windows as well as sharp fluctuations in GPP and LE, respectively. This error in prediction of fluxes ($\mathrm{\Delta}=\mathrm{Observed}\mathrm{Posterior}$) is also revealed by posterior simulations represented by brown and orange lines. Comparing the posterior simulations with and without the MODIS reflectance constraints with the prior reveals the net error to be of similar magnitude in both cases with δ_{Δ}≤15 % during the growing season and indicates equifinality in the predicted fluxes with the parameter combinations.
Figure 11 shows the final posterior error reduction of the retrieval iterations for each moving window. It is found that there is significant reduction in the posterior errors for BB_{slope}. We find that ${\mathit{\zeta}}_{{V}_{\mathrm{cmax}}}$ and ζ_{LAI} have the same trend and seasonality as that of retrieved V_{cmax} and LAI respectively. The evolution of the retrieval error correlations is again interesting and comparison with the Mead CornC_{4} site using flux observations only shows it is similar for ${\mathit{\rho}}_{{V}_{\mathrm{cmax},\mathrm{LAI}}}$ and ${\mathit{\rho}}_{{\mathrm{BB}}_{\mathrm{slope},\mathrm{LAI}}}$ (negative and positive respectively) and different for ${\mathit{\rho}}_{{V}_{\mathrm{cmax}},{\mathrm{BB}}_{\mathrm{slope}}}$.
The correlations ${\mathit{\rho}}_{{V}_{\mathrm{cmax}},{\mathrm{BB}}_{\mathrm{slope}}}$ and ${\mathit{\rho}}_{{V}_{\mathrm{cmax},\mathrm{LAI}}}$ are both negative, indicating counteracting effects of these variables in the retrieval. In comparison, ${\mathit{\rho}}_{{\mathrm{BB}}_{\mathrm{slope},\mathrm{LAI}}}$ is positive for the retrievals, indicating insync behavior. The error correlations ${\mathit{\rho}}_{{V}_{\mathrm{cmax},\mathrm{LAI}}}$ are highest among the three in the middle of the growing season.
Finally, there is a significant improvement in the estimation of both GPP and LE fluxes (see Fig. S6) (R^{2}=0.7) with the optimized state vector over the prior values. The optimized state vector is able to capture the growing seasonal variability well (slope of regression lines: GPP =0.97 and LE =1.02). The unoptimized prior values severely underpredict both fluxes. The corresponding error measures ${\mathit{\chi}}_{\mathrm{GPP}\mathrm{opt}}^{\mathrm{2}}=\mathrm{12}\phantom{\rule{0.125em}{0ex}}\mathrm{975}$, ${\mathit{\chi}}_{\mathrm{GPP}\mathrm{unopt}}^{\mathrm{2}}=\mathrm{30}\phantom{\rule{0.125em}{0ex}}\mathrm{070}$, ${\mathit{\chi}}_{\mathrm{LE}\mathrm{opt}}^{\mathrm{2}}=\mathrm{17}\phantom{\rule{0.125em}{0ex}}\mathrm{260}$ and ${\mathit{\chi}}_{\mathrm{LE}\mathrm{unopt}}^{\mathrm{2}}=\mathrm{32}\phantom{\rule{0.125em}{0ex}}\mathrm{283}$ indicate that the inversion framework is able to retrieve the seasonal dynamics in the photosynthetic and canopy structural parameters for accurate prediction of the fluxes.
6.4.3 Inversion results for the year 2007
The year 2007 for the Missouri Ozark site is an interesting example because the forest experiences a late summer drought and decrease in productivity, which is captured by the flux observations but not with phenology or greenness (see Fig. S5, Figs. 12 and 5). The inversion setup, meteorological data resolution, initial guess and prior errors for the year 2007 are similar to that of 2009. Figure 12 shows the results for the retrieval of parameters V_{cmax}, BB_{slope} and LAI from the inversion framework using (i) observations of GPP and LE fluxes (dashed line) and (ii) flux observations with MODIS reflectance bands (solid lines). The overall seasonality of the retrieved parameters reveals that the inversion framework is able to capture the latesummer drought. The midseason drought around DOY 230–250 is well captured by decreases in photosynthesis (V_{cmax}) and stomatal conductance (BB_{slope}). These parameters again increase from around DOY 260 after the drought recovery and also correlate well with the flux observations. However, the V_{cmax} variations should not be confused with actual changes in rubisco concentrations. Like many other carbon cycle models, the only way to impose environmental stress (apart from VPD and temperature) in SCOPE is to scale V_{cmax} with a stress factor between 0 and 1. Our retrieved effective V_{cmax} is the product of a stress factor and the physiological V_{cmax}, which explains the large variations during drought here.
In contrast, the change in retrieved LAI during the period from DOY 240 to 300 is fairly gradual and reflects the phenology only. The addition of MODIS reflectance data constrains LAI (and in conjunction V_{cmax}) much better than just the flux observations, which is clearly evident from the V_{cmax} retrievals around DOY 200. This large change in V_{cmax} (when using just constraining flux observations) also corresponds to the single largest rainfall event of the season (see Fig. S5) as well as a corresponding concurrent spike in productivity. The MODIS red and NIR reflectances, unlike GPP and LE fluxes, were insensitive to this short drought and provide better constraints on the LAI and V_{cmax} retrievals during this period. A comparison of modeled red and NIR reflectance from the SCOPE model with optimized parameters shows that it matches well with the observations (see Fig. S8). There is excellent match in the red reflectance throughout the season, and the NIR reflectance also matches well with the observations during the early–middle part of the growing season (DOYs 130–250); during the postdrought recovery phase the increase may be attributed to the increase in retrieved LAI. Apart from the drought period, the range and variability of all three parameters correspond well with the retrievals for the same site for the year 2009 presented earlier.
The posterior simulations excluding and including the MODIS observations (represented as orange and brown lines respectively in Fig. 12) indicate similar improvement in prediction of fluxes over their priors with ${R}_{\mathrm{GPP}\mathrm{opt}}^{\mathrm{2}}=\mathrm{0.7}$, ${R}_{\mathrm{GPP}\mathrm{unopt}}^{\mathrm{2}}=\mathrm{0.2}$, ${R}_{\mathrm{LE}\mathrm{opt}}^{\mathrm{2}}=\mathrm{0.5}$ and ${R}_{\mathrm{LE}\mathrm{unopt}}^{\mathrm{2}}=\mathrm{0.1}$. For the year 2007 it is also found that the posterior error reduction ${\mathit{\zeta}}_{{V}_{\mathrm{cmax}}}$ and ζ_{LAI} is similar or slightly better compared to the year 2009 (see Fig. S7). This again provides more confidence in the retrieved parameters and their temporal dynamics. As expected, the posterior error reduction is slightly better for all three parameters with the reflectance constraints during the middle of the growing season. The error correlations ${\mathit{\rho}}_{{V}_{\mathrm{cmax}},{\mathrm{BB}}_{\mathrm{slope}}}$ and ${\mathit{\rho}}_{{V}_{\mathrm{cmax},\mathrm{LAI}}}$ for 2007 are both negative and similar in magnitude to that of year 2009. The error correlation structure for the year 2007 is different compared to the wet year 2009, especially during the middle of the growing season, and may be attributed to the change in interaction between the state vector elements due to imposed drought stresses. We note that our moving window inversion setup using different observational streams is able to capture the ecosystem dynamics over the entire growing season as well as capture inseason drought dynamics, which are not possible using traditional onestep seasonal or annual inversion approaches. Inversions like this will also help guide model parameterizations of stress impacts on the dynamic downregulation of photosynthesis as a response to, for example, changes in the soil matric potential.
Finally, we performed sensitivity analyses of the newer temperature dependence implementation in SCOPE on the inversion retrieval results for the Ozark site (see Sect. S1.1.2). It is found that with the newer temperature dependence implementation and optimized parameters SCOPE is able to capture V_{cmax} variations due to changes in the average canopy temperatures well. There is a clear difference between the retrieved V_{cmax} with and without temperature dependence, and the changes correlate with the implemented temperature response functions (see Sect. S1.1.2). The results of the posterior simulation on GPP and LE fluxes also indicate improvement with δχ^{2} error reduction of 3415 and 2104 for GPP and LE respectively.
6.5 Retrieval results for the Niwot Ridge site
6.5.1 Site description
The Niwot Ridge site is also a part of the AmeriFlux network located in a subalpine forest ecosystem just below the continental divide near Nederland, Colorado. The average elevation of this site is 3050 m and it is one of the high alpine evergreen needleleaf forests with C_{3} plant species (Burns et al., 2016). This ecosystem is nearly 100 years old, and thus very different from the Mead and the Ozark sites (Monson et al., 2002). This site has a mean annual precipitation of 800 mm and a mean annual temperature of 1.5 ^{∘}C and has a continuous data record from 1998 until the present. This site is thus the coldest and driest among the three. We choose the year 2010 for the current analysis, using halfhourly flux data. The site meteorology and forcing variables relevant to the SCOPE inversion retrievals for the year 2010 are presented in the Fig. S11. The snow cover at this site affects the energy exchanges and GPP for a large fraction of the year. The trees are evergreen and there is no well defined growing season but we find that most photosynthetic activity occurs in the period between May and October, with a smaller flux magnitude compared to either the Mead or Ozark sites. The sensible heat at the site is also larger during this period compared to the latent heat fluxes. Since the LAI for the site does not really change over the year, we focus on the retrieval of the parameters V_{cmax} and BB_{slope} only from constraining GPP and LE fluxes, fixing the LAI.
6.5.2 Inversion parameters and results
For the retrievals the prior value of the state vector along with prior errors and daytime duration used in the retrieval windows are shown in Table 1. For this evergreen site, we assumed a prior value of LAI equal to 3.8 (Monson et al., 2009) with very low variance, effectively fixing the LAI, which improves the retrieval of the other state vector parameters as it reduces error covariations.
Figure 13 shows the results for the retrieval of parameters V_{cmax}, BB_{slope} and LAI. The grey time series of GPP and LE in the background are the actual values used for constructing the observation vector Y corresponding to each retrieval window for parameter retrieval. The V_{cmax} values again follow similar trends to the GPP fluxes across the entire active season and the inversion captures the rise and fall in the GPP trends extremely well. The slight dip and rise in the DOY 190–210 period follows the GPP observational data and may be attributed to temperature and precipitation fluctuations. The BB_{slope} seems to closely follow variations in LE fluxes and captures the seasonality well. However, we should note that changes in soil evaporation might alias into the retrieval of the BB_{slope}, especially at Niwot Ridge.
Figure 14 shows the final posterior error reduction of the retrieval iterations for each moving window. The posterior error reductions for BB_{slope} and V_{cmax} are high throughout the year, likely attributable to a constant LAI used for this example (e.g. at vanishing LAI in winter for deciduous ecosystems, the Jacobians of fluxes with respect to physiological parameters will be meaningless and vanish as well, thereby reducing error reductions). The evolution of the error correlations follows similar trends to those in the Missouri Ozark C_{3} site. We find ${\mathit{\rho}}_{{V}_{\mathrm{cmax}},{\mathrm{BB}}_{\mathrm{slope}}}$ and ${\mathit{\rho}}_{{V}_{\mathrm{cmax},\mathrm{LAI}}}$ are both negative and ${\mathit{\rho}}_{{\mathrm{BB}}_{\mathrm{slope},\mathrm{LAI}}}$ is positive. There is a sharp discontinuity around DOY 250–260 in terms of ${\mathit{\zeta}}_{{\mathrm{BB}}_{\mathrm{slope}}}$ and ${\mathit{\rho}}_{{V}_{\mathrm{cmax},\mathrm{LAI}}}$, ${\mathit{\rho}}_{{V}_{\mathrm{cmax}},{\mathrm{BB}}_{\mathrm{slope}}}$, probably due to quality and/or discontinuity in the observational fluxes and environmental forcings. The prior and posterior simulations (with unoptimized and optimized parameters respectively) show overall net improvement of the flux predictions when compared to tower observations with ${R}_{\mathrm{GPP}\mathrm{opt}}^{\mathrm{2}}=\mathrm{0.85}$, ${R}_{\mathrm{GPP}\mathrm{unopt}}^{\mathrm{2}}=\mathrm{0.79}$, ${R}_{\mathrm{LE}\mathrm{opt}}^{\mathrm{2}}=\mathrm{0.52}$ and ${R}_{\mathrm{LE}\mathrm{unopt}}^{\mathrm{2}}=\mathrm{0.49}$. The corresponding error measures are ${\mathit{\chi}}_{\mathrm{GPP}\mathrm{opt}}^{\mathrm{2}}=\mathrm{38}\phantom{\rule{0.125em}{0ex}}\mathrm{971}$, ${\mathit{\chi}}_{\mathrm{GPP}\mathrm{unopt}}^{\mathrm{2}}=\mathrm{22}\phantom{\rule{0.125em}{0ex}}\mathrm{085}$, ${\mathit{\chi}}_{\mathrm{LE}\mathrm{opt}}^{\mathrm{2}}=\mathrm{63}\phantom{\rule{0.125em}{0ex}}\mathrm{865}$ and ${\mathit{\chi}}_{\mathrm{LE}\mathrm{unopt}}^{\mathrm{2}}=\mathrm{68}\phantom{\rule{0.125em}{0ex}}\mathrm{463}$. The high χ^{2} values for posterior GPP over the prior may indicate that the posterior simulation may not capture some periods in the simulation where the observed fluxes are higher than average (which were probably captured with a high prior value of V_{cmax}). This may also be due to model structural issues caused by stress induced due to snow processes not being accurately represented in SCOPE (with a stress factor in V_{cmax}). Further, we have an additional constraint of nearconstant LAI, while we have found that GPP is very sensitive to LAI. In this case the inversion has to fit GPP by mainly varying the V_{cmax}, which gives it fewer degrees of freedom and may not yield the most optimal solution. For instance, changes in the fraction of direct vs. diffuse light is not fully represented in our model (apart from changing PAR) but could affect the overall PAR value incident at the top of the canopy as well as the overall light interception, as there are considerable gaps between canopies.
Our results demonstrate the feasibility of a moving window inversion approach for the retrieval of key ecosystem parameters using eddy covariance flux tower observations. In addition, we also demonstrated that red and NIR spectral reflectance observation from satellites adds better constraints on LAI, and thereby also improves the retrievals of V_{cmax} and BB_{slope} by reducing interferences in retrieved parameters. The moving window retrieval approach is specifically useful for dynamic changes in ecosystem parameters, such as the response to environmental stress due to water stress, which we observed during a summer drought at the Ozark flux tower site.
There is strong evidence from measurements that under normal conditions both LAI and photosynthetic parameters have seasonal variability (Wang et al., 2008; Wilson and Baldocchi, 2000; Wilson et al., 2000), which correlate with energy fluxes. Our model inversion results are in alignment and agree well with these findings. From our results, we find considerable seasonal variability in V_{cmax} and BB_{slope} (to some extent). Previously, many studies have reported measured V_{cmax} values of similar ranges, such as 0–70 µmols m^{−2} s^{−1} (Wilson et al., 2000) for deciduous trees and 0–80 µmols m^{−2} s^{−1} annual variability in tall Japanese red pine forests (Han et al., 2004). In addition, most of the V_{cmax} variability is found in systems which have seasonally variable or constant nitrogen (N) content (Wilson et al., 2000). These changes may be mostly attributed to substantial inseason changes in the fraction of total N allocated to rubisco as well as changes in leaf mass per area (Wilson et al., 2000). In addition, most models, including SCOPE, have no other method of imposing environmental stress than reducing V_{cmax} by a stress factor [0, 1]. The effect of reductions in V_{cmax} are a reduction in assimilation, which also suppresses transpiration. Thus, we are fitting an effective V_{cmax} parameter, which factors in effects from true changes in rubisco content as well as the impact of environmental stress. It has been shown that there is also a possibility to have large seasonality variability in BB_{slope} (15 to 25) values (Wolf et al., 2006). The SCOPE model handles both the C_{3} and C_{4} photosynthetic pathways and could thus be applied to study a wide variety of ecosystems. We demonstrate the approach here for climate and productivity gradients across agricultural, deciduous broadleaf forest and subalpine evergreen forest ecosystems.
The framework demonstrates the feasibility of the approach in parameter retrievals using suitable a priori uncertainty in state vector and observation noise. The uncertainty in surface energy balance closure has been generally reported to be around 10 %–30 % (Sánchez et al., 2010; Wilson et al., 2002; von Randow et al., 2004) and is found to be dependent on timescales due to differences in energy storage terms in ecosystems (Reed et al., 2018). Apart from observations, it should also be noted that filtering and quality control may be necessary for the meteorological and/or forcing fluxes, as artifacts in the data can influence the inversion and optimization and greatly affect the results.
In the current implementation, the inversion approach may not properly retrieve the key parameters for ecosystems for which the underlying forward model (SCOPE) may have deficiencies in process representation. For example, there are competing optimality theories between whether the BB_{slope} (Mäkelä et al., 1996; Van der Tol et al., 2008a, b) or V_{cmax} (Xu and Baldocchi, 2003) is most affected by drought during the growing season. An improvement in the current framework in the form of better process representation (and parameterization) of the soil moisture status in the stomatal conductance model within SCOPE may improve the results and could serve as a test bed for improvements in process representations. In this case, it could be through the implementation of leaf water potential (Tuzet et al., 2003) or an optimality approach between water loss and carbon gain (Katul et al., 2009; Medlyn et al., 2011).
While this paper focuses on the conceptual inversion framework, demonstrating a novel approach for estimating important ecosystem parameters in short time windows for modeling the dynamics of coupled carbon and water fluxes across ecosystems, there are opportunities for improving the overall inversion approach to better estimate the parameters. SCOPE allows us to ingest a variety of other observables to constrain the parameter space, including spectrally resolved reflectance (which can constrain LAI and chlorophyll content) as well as thermal emissions (which constrain LE) and SIF (which constrains APAR and V_{cmax}). Global time series of SIF observations, which provide a direct probe into photosynthetic machinery, are now available from spacebased (Frankenberg et al., 2011, 2014; Guanter et al., 2014) and groundbased observations (Frankenberg et al., 2016). It will also be important to quantify the respective information content for the observables within the framework, which can vary depending on the state vector itself as the Jacobians from our inversion results indicate that the inverse system is nonlinear.
Our Bayesian inversion framework is highly flexible in terms of allowing an arbitrary number of prior and retrieval parameters, which could be tuned for better estimation of the key ecosystem parameters with accurate posterior uncertainty estimates. The stepwise LM optimization framework with the SCOPE model also automatically weighs the carbon and water fluxes towards optimal state vector estimation without any predefined constraints (Wolf et al., 2006) for particular parameters. Further, the Bayesian retrieval framework could also be used to retrieve the photosynthetic temperature dependency parameters such as entropies and activation energies, which are even harder to measure directly but might be crucial, especially as the modeling of the ecosystem response to a warming climate mostly neglects potential changes in temperature dependencies of V_{cmax}. Our ongoing and future research efforts aim to utilize the developed framework to address these questions by incorporating newer observations and making use of the full potential of SCOPE.
The authors thank the AmeriFlux team for making the eddycovariance flux data available for this study. The FLUXNET2015 datasets used in this study have been downloaded from the FLUXNET community data portal (http://fluxnet.fluxdata.org/data/fluxnet2015dataset/, FLUXNET community, 2015). The version of SCOPE model used in this study can be obtained from https://github.com/Christiaanvandertol/SCOPE (van der Tol, 2017).
The biochemical module is at the center of energy balance computations within SCOPE. This module computes the net assimilation (photosynthesis), stomatal conductance and the chlorophyll fluorescence of a leaf. This module is thus extremely important because the coupled photosynthesis and stomatal conductance regulates the latent heat flux which in turn affects the net energy balance and the leaf temperature, which in turn again affects the leaf photosynthesis and subsequently the energy balance. SCOPE computes the leaf temperature and the overall energy balance iteratively such that they there is closure in energy balance. As such the leaf temperature and its regulation of photosynthesis forms an extremely important component of the overall energy balance of the canopy. We have adapted a photosynthetic model together with coupled temperature dependence of the photosynthetic parameters according to the implementation in CLM4.5. This includes both the temperature dependence functions and the high temperature inhibition of the parameters. The model includes exclusive pathways both for the C_{3} and C_{4} plant species and is represented as follows:
The net photosynthesis (assimilation) after accounting for respiration (R_{d}) is given as follows:
Further, the rate limiting steps are represented as follows. The RuBP carboxylase (rubisco) limited rate of carboxylation A_{c} is given by they following:
The lightlimited rate of carboxylation (governed by the capacity to regenerate RuBP) A_{j} is given by the following:
Finally the productlimited carboxylation rate for C_{3} plants and the PEPcarboxylaselimited rate of carboxylation for the C_{4} plants A_{p} is given by the following:
For the above Eqs. A2, A3 and A4, we have the assimilation rates ${A}_{\mathrm{c},j,\mathrm{p}}$ (in the units of µmols m^{−2} s^{−1}), C_{i} is the internal CO_{2} concentration of the leaf (units of ppm) and V_{cmax} is the maximum rate of carboxylation. For the C_{3} species, K_{c} and K_{o} are the Michaelis–Menten constants for CO_{2} and O_{2} respectively (units of µmols m^{−2} s^{−1}), Γ^{*} is the CO_{2} compensation point (units are ppm), J is the potential electron transport rate (units of µmols m^{−2} s^{−1}) and T_{p} is the triose phosphate utilization rate. For the C_{4} plants, ϕ is the absorbed PAR (in the units of Wm^{−2}) and the factor 4.6 converts it to PPFD (photosynthetic photon flux density, in units of µmol m^{−2} s^{−1}) (for SCOPE biochem module the PAR is already in PPFD units), α is the quantum efficiency (0.05 mol CO_{2} mol^{−1} photon), and k_{p} is the initial slope of the C_{4} CO_{2} response curve.
For the C_{3} plants, the potential electron transport rate J depends on the PAR absorbed by a leaf, which is obtained as the smaller root of the two roots of the following equation:
where J_{max} is the maximum electron transport rate (µmols m^{−2} s^{−1}), I_{PSII} is the light used in photosystem II (µmols m^{−2} s^{−1}) which is given by Eq. (A6) and Θ_{PSII} is a curvature parameter.
The term Φ_{PSII} in Eq. (A6) is the quantum yield of photosystem II and 0.5 represents halfelectron transfer to each of the photosystems I and II. The overall gross photosynthesis rate is computed as a colimitation (Collatz et al., 1991a, 1992) and is computed as the smaller root of the following equations:
The parameters Θ_{cj} and Θ_{ip} control the smoothness of the light response curve between lightlimited and enzyme or productlimiting rates. The values of the different parameters at optimum temperature (mostly as a function of V_{cmax25}, here the optimum temperature is assumed to be 25 ^{∘}C) used in the photosynthesis model are presented in Table A1.
The photosynthesis model parameters for both the C_{3} and C_{4} pathways described in the previous section and shown in Table A1 have temperaturedependent variations and need to be adjusted for specific leaf temperature before implementing them in the photosynthesis model. The temperature dependence of photosynthetic parameters for the C_{3} species can be broadly decomposed into two parts: (i) the temperature response and (ii) the high temperature inhibition. The functional form of these are as follows:
where H_{a} is the activation energy, H_{d} is the deactivation energy, S_{v} is the entropy term, T_{0} is the optimum temperature and T_{v} is the leaf temperature. The functional relationship of the different photosynthetic parameters in the C_{3} pathway are as follows:
The temperature dependence of photosynthetic parameters for the C_{4} species are given by the following relationships:
The Q_{10} temperature coefficient is a measure of the rate of change of a biological or chemical system as a consequence of increasing the temperature by 10 ^{∘}C. The values of the temperature dependence functional parameters for both C_{3} and C_{4} species used in the present study are provided in Tables B1 and B2 respectively.
The temperature dependence parameters (activation, deactivation and entropy) are variable between different plant species (Leuning, 2002) as such its formulation in the newer implementation of the SCOPE model allows us to use appropriate values depending on the ecosystem we study.
Ag_{s}C_{i} iterations
The final solution for photosynthesis requires an iterative solution of the coupled equations representing (i) the Farquhar, von Caemmerer and Berry (FvCB) model (Farquhar et al., 1980) for the photosynthesis rate (A); (ii) Fick's law of diffusion (Eq. B6) for internal (C_{i}) and leaf surface (C_{s}) CO_{2} concentration; and (iii) Ball–Berry stomatal conductance model (Ball et al., 1987) (Eq. B7) for stomatal conductance (g_{s}) to obtain stable converging solutions.
In Eq. (B7), BB_{slope} represents the Ball–Berry slope, r_{h} the relative humidity and g_{0} the Ball–Berry intercept. In the absence of an initial specification of C_{i}, we make the assumption that g_{0}=0 in Eq. (B7); then, combining Eqs. (B6), (B7), the initial estimate of C_{i} is given as follows:
where ${f}_{C\mathrm{i}}^{\mathrm{min}}$ is the assumed minimum fractional leaf boundary layer CO_{2} (assumed to be 0.3 for C_{3} and 0.1 for C_{4} species). This initial estimate of C_{i} is used to again estimate the photosynthesis based on the FvCB model (Farquhar et al., 1980), followed by estimation of stomatal conductance using the Ball–Berry model Eq. (B7). Finally, the Newton–Raphson method is used to obtain a forward estimation of the new value of internal CO_{2} concentration (Ivanov et al., 2008; Sun et al., 2012). The updated C_{i} is further used in the Ag_{s}C_{i} until convergence.
For deriving the maximum probability state X ($\widehat{\mathit{X}}$) we equate the derivative of the Eq. (5) to zero to obtain the following:
It can be noted here that the gradient ∇_{X} of the above vectorvalued function is a matrixvalued function and the Jacobian matrix is represented as K(X)=∇_{X}F(X), which results in the following implicit equation for $\widehat{\mathit{X}}$:
We have to now use any general rootfinding method for finding the solutions of Eq. (C2). If the problem is not too nonlinear we can use the Newton and Gauss–Newton iterative methods (Hartley, 1961). In general for any vector equation G(X)=0, we can write the Newton iteration as follows:
For our problem we can assume the derivative of the costfunction G(X) to be the lefthand side of Eq. (C1); therefore the gradient of G(X) (∇G) also known as the Hessian is given by the following:
The Hessian in Eq. (C4) involves the Jacobian K and both the first and second derivatives of the forward model. The second derivative is complicated because it is a vector whose elements are matrices and further this term is postmultiplied by the factor ${\mathbf{S}}_{\mathit{\epsilon}}^{\mathrm{1}}[\mathit{Y}F(\mathit{X}\left)\right]$. The third term in the righthand side of Eq. (C4) is thus computationally expensive and further for moderately linear problems this term is small; as such this term can be ignored (also called smallresidual problems in numerical methods). When we ignore this term, we get the Gauss–Newton iteration scheme by substituting Eqs. (C2) and (C4) in Eq. (C3):
where K_{i}=K(X_{i}). We can substitute F(X) from Eq. (4) in Eq. (C2) to get the following:
Again, representing ${\mathbf{K}}_{\mathrm{l}}={\mathrm{\nabla}}_{X}F(\mathit{X}{)}_{\mathit{X}={\mathit{X}}_{\mathrm{l}}}$, ${\mathit{F}}_{\mathrm{l}}=F(\mathit{X}{)}_{\mathit{X}={\mathit{X}}_{\mathrm{l}}}$ we can further simplify and rearrange Eq. (C6) as follows:
In the above equations, if we change the interpretation of the subscript l from “linearization” to “iteration counter”, we obtain the following equation:
If we express X_{i+1} as a departure from X_{i} rather than X_{a} we obtain the same expression for the iteration steps as Eqs. (C5) or (6).
The supplement related to this article is available online at: https://doi.org/10.5194/bg16772019supplement.
DD and CF designed the research, developed the inversion framework and conducted the modeling simulations. DD, CF, CvdT, YS and DSS contributed towards interpretation and analysis of the results. DD and CF developed the initial draft of the manuscript, and DD, CF, CvdT, YS and DSS contributed towards improvement of the original draft and finalization of the paper.
The authors declare that they have no conflict of interest.
AmeriFlux site Missouri Ozark (USMOz) is supported by the US Department of Energy (DOE), Office of Science, Office of Biological and Environmental Research Program, through Oak Ridge National Laboratory's Terrestrial Ecosystem ScienceScience Focus Area; ORNL is managed by UTBattelle, LLC, for DOE under contract DEAC0500OR22725.
The research was carried out, in part, at the Jet Propulsion Laboratory,
California Institute of Technology, under a contract with the National
Aeronautics and Space Administration.
Government sponsorship is acknowledged.
Edited by: Dan Yakir
Reviewed by: Peter Rayner and two
anonymous referees
Baldocchi, D.: An analytical solution for coupled leaf photosynthesis and stomatal conductance models, Tree Physiol., 14, 1069–1079, 1994. a
Baldocchi, D., Valentini, R., Running, S., Oechel, W., and Dahlman, R.: Strategies for measuring and modelling carbon dioxide and water vapour fluxes over terrestrial ecosystems, Glob. Change Biol., 2, 159–168, 1996. a
Baldocchi, D., Kelliher, F. M., Black, T. A., and Jarvis, P.: Climate and vegetation controls on boreal zone energy exchange, Glob. Change Biol., 6, 69–83, 2000. a
Baldocchi, D., Falge, E., Gu, L., Olson, R., Hollinger, D., Running, S., Anthoni, P., Bernhofer, C., Davis, K., Evans, R., et al.: FLUXNET: A new tool to study the temporal and spatial variability of ecosystem–scale carbon dioxide, water vapor, and energy flux densities, B. Am. Meteorol. Soc., 82, 2415–2434, 2001. a
Ball, J. T., Woodrow, I. E., and Berry, J. A.: A model predicting stomatal conductance and its contribution to the control of photosynthesis under different environmental conditions, Progress in photosynthesis research, 221–224, Springer, Dordrecht, 1987. a, b
Beerling, D. and Quick, W.: A new technique for estimating rates of carboxylation and electron transport in leaves of C3 plants for use in dynamic global vegetation models, Glob. Change Biol., 1, 289–294, 1995. a
Bernacchi, C. J., Singsaas, E. L., Pimentel, C., Portis Jr., A. R., and Long, S. P.: Improved temperature response functions for models of Rubiscolimited photosynthesis, Plant Cell Environ., 24, 253–259, 2001. a, b
Bonan, G. B., Lawrence, P. J., Oleson, K. W., Levis, S., Jung, M., Reichstein, M., Lawrence, D. M., and Swenson, S. C.: Improving canopy processes in the Community Land Model version 4 (CLM4) using global flux fields empirically inferred from FLUXNET data, J. Geophys. Res.Biogeo., 116, https://doi.org/10.1029/2010JG001593, 2011. a, b, c, d
Burns, S. P., Maclean, G. D., Blanken, P. D., Oncley, S. P., Semmer, S. R., and Monson, R. K.: The Niwot Ridge Subalpine Forest USNR1 AmeriFlux site – Part 1: Data acquisition and site recordkeeping, Geosci. Instrum. Method. Data Syst., 5, 451–471, https://doi.org/10.5194/gi54512016, 2016. a
Chen, J. M., Rich, P. M., Gower, S. T., Norman, J. M., and Plummer, S.: Leaf area index of boreal forests: Theory, techniques, and measurements, J. Geophys. Res.Atmos., 102, 29429–29443, 1997. a, b
Collatz, G., Ball, J., Grivet, C., and Berry, J. A.: Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: a model that includes a laminar boundary layer, Agr. Forest Meteorol., 54, 107–136, https://doi.org/10.1016/01681923(91)900028, 1991a. a
Collatz, G. J., Ball, J. T., Grivet, C., and Berry, J. A.: Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: a model that includes a laminar boundary layer, Agr. Forest Meteorol., 54, 107–136, 1991b. a
Collatz, G. J., RibasCarbo, M., and Berry, J.: Coupled photosynthesisstomatal conductance model for leaves of C4 plants, Funct. Plant Biol., 19, 519–538, 1992. a, b
Cox, P. M., Pearson, D., Booth, B. B., Friedlingstein, P., Huntingford, C., Jones, C. D., and Luke, C. M.: Sensitivity of tropical carbon to climate change constrained by carbon dioxide variability, Nature, 494, p. 341, 2013. a
Cramer, W., Bondeau, A., Woodward, F. I., Prentice, I. C., Betts, R. A., Brovkin, V., Cox, P. M., Fisher, V., Foley, J. A., Friend, A. D., et al.: Global response of terrestrial ecosystem structure and function to CO_{2} and climate change: results from six dynamic global vegetation models, Glob. Change Biol., 7, 357–373, 2001. a
Dai, Y., Dickinson, R. E., and Wang, Y.P.: A twobigleaf model for canopy temperature, photosynthesis, and stomatal conductance, J. Climate, 17, 2281–2299, 2004. a
Dutta, D., Wang, K., Lee, E., Goodwell, A., Woo, D. K., Wagner, D., and Kumar, P.: Characterizing Vegetation Canopy Structure Using Airborne Remote Sensing Data, IEEE T. Geosci. Remote, 55, 1160–1178, 2017. a, b
Falkowski, P., Scholes, R. J., Boyle, E. E. A., Canadell, J., Canfield, D., Elser, J., Gruber, N., Hibbard, K., Högberg, P., Linder, S., and Mackenzie, F. T.: The global carbon cycle: a test of our knowledge of earth as a system, Science, 290, 291–296, 2000. a
Farquhar, G. v., von Caemmerer, S. V., and Berry, J.: A biochemical model of photosynthetic CO_{2} assimilation in leaves of C3 species, Planta, 149, 78–90, 1980. a, b, c, d, e
Flexas, J., Escalona, J. M., Evain, S., Gulías, J., Moya, I., Osmond, C. B., and Medrano, H.: Steadystate chlorophyll fluorescence (Fs) measurements as a tool to follow variations of net CO_{2} assimilation and stomatal conductance during waterstress in C3 plants, Physiol. Plantarum, 114, 231–240, 2002. a
FLUXNET community: FLUXNET 2015 dataset, available at: http://fluxnet.fluxdata.org/data/fluxnet2015dataset/, last access: 26 October 2015.
Frankenberg, C., Fisher, J. B., Worden, J., Badgley, G., Saatchi, S. S., Lee, J. E., Toon, G. C., Butz, A., Jung, M., Kuze, A., and Yokota, T.: New global observations of the terrestrial carbon cycle from GOSAT: Patterns of plant fluorescence with gross primary productivity, Geophys. Res. Lett., 38, https://doi.org/10.1029/2011GL048738, 2011. a, b
Frankenberg, C., O'Dell, C., Berry, J., Guanter, L., Joiner, J., Köhler, P., Pollock, R., and Taylor, T. E.: Prospects for chlorophyll fluorescence remote sensing from the Orbiting Carbon Observatory2, Remote Sens. Environ., 147, 1–12, 2014. a
Frankenberg, C., Drewry, D., Geier, S., Verma, M., Lawson, P., Stutz, J., and Grossmann, K.: Remote sensing of solar induced Chlorophyll Fluorescence from satellites, airplanes and groundbased stations, in: Geoscience and Remote Sensing Symposium (IGARSS), 2016 IEEE International, 1707–1710, IEEE, 2016. a
Friedlingstein, P., Cox, P., Betts, R., Bopp, L., von Bloh, W., Brovkin, V., Cadule, P., Doney, S., Eby, M., Fung, I., and Bala, G.: Climate–carbon cycle feedback analysis: results from the C4MIP model intercomparison, J. Climate, 19, 3337–3353, 2006. a
Gu, L., Meyers, T., Pallardy, S. G., Hanson, P. J., Yang, B., Heuer, M., Hosman, K. P., Riggs, J. S., Sluss, D., and Wullschleger, S. D.: Direct and indirect effects of atmospheric conditions and soil moisture on surface energy partitioning revealed by a prolonged drought at a temperate forest site, J. Geophys. Res.Atmos., 111, https://doi.org/10.1029/2006JD007161, 2006. a
Guanter, L., Zhang, Y., Jung, M., Joiner, J., Voigt, M., Berry, J. A., Frankenberg, C., Huete, A. R., ZarcoTejada, P., Lee, J. E., and Moran, M. S.: Global and timeresolved monitoring of crop photosynthesis with chlorophyll fluorescence, P. Natl. Acad. Sci. USA, 111, E1327–E1333, 2014. a
Han, Q., Kawasaki, T., Nakano, T., and Chiba, Y.: Spatial and seasonal variability of temperature responses of biochemical photosynthesis parameters and leaf nitrogen content within a Pinus densiflora crown, Tree Physiol., 24, 737–744, 2004. a
Hartley, H. O.: The modified GaussNewton method for the fitting of nonlinear regression functions by least squares, Technometrics, 3, 269–280, 1961. a
Heimann, M. and Reichstein, M.: Terrestrial ecosystem carbon dynamics and climate feedbacks, Nature, 451, 7176, https://doi.org/10.1038/nature06591, 2008. a
Houborg, R., Soegaard, H., and Boegh, E.: Combining vegetation index and model inversion methods for the extraction of key vegetation biophysical parameters using Terra and Aqua MODIS reflectance data, Remote Sens. Environ., 106, 39–58, 2007. a
Ivanov, V. Y., Bras, R. L., and Vivoni, E. R.: Vegetationhydrology dynamics in complex terrain of semiarid areas: 1. A mechanistic approach to modeling dynamic feedbacks, Water Resour. Res., 44, https://doi.org/10.1029/2006WR005588, 2008. a
Jacquemoud, S., Baret, F., Andrieu, B., Danson, F., and Jaggard, K.: Extraction of vegetation biophysical parameters by inversion of the PROSPECT+ SAIL models on sugar beet canopy reflectance data. Application to TM and AVIRIS sensors, Remote Sens. Environ., 52, 163–172, 1995. a
Katul, G., Manzoni, S., Palmroth, S., and Oren, R.: A stomatal optimization theory to describe the effects of atmospheric CO_{2} on leaf photosynthesis and transpiration, Ann. Bot.London, 105, 431–442, 2009. a
Knorr, W. and Heimann, M.: Uncertainties in global terrestrial biosphere modeling: 1. A comprehensive sensitivity analysis with a new photosynthesis and energy balance scheme, Global Biogeochem. Cy., 15, 207–225, 2001. a
Knorr, W. and Kattge, J.: Inversion of terrestrial ecosystem model parameter values against eddy covariance measurements by Monte Carlo sampling, Glob. Change Biol., 11, 1333–1351, 2005. a
Kucharik, C. J., Foley, J. A., Delire, C., Fisher, V. A., Coe, M. T., Lenters, J. D., YoungMolling, C., Ramankutty, N., Norman, J. M., and Gower, S. T.: Testing the performance of a dynamic global ecosystem model: water balance, carbon balance, and vegetation structure, Global Biogeochem. Cy., 14, 795–825, 2000. a
Lawrence, D. M., Oleson, K. W., Flanner, M. G., Thornton, P. E., Swenson, S. C., Lawrence, P. J., Zeng, X., Yang, Z. L., Levis, S., Sakaguchi, K., and Bonan, G. B.: Parameterization improvements and functional and structural advances in version 4 of the Community Land Model, J. Adv. Model. Earth Sy., 3, https://doi.org/10.1029/2011MS00045, 2011. a
Leuning, R.: Temperature dependence of two parameters in a photosynthesis model, Plant Cell Environ., 25, 1205–1210, 2002. a, b, c, d, e
Leuning, R., Van Gorsel, E., Massman, W. J., and Isaac, P. R.: Reflections on the surface energy imbalance problem, Agr. Forest Meteorol., 156, 65–74, 2012. a
Levenberg, K.: A method for the solution of certain nonlinear problems in least squares, Q. Appl. Math., 2, 164–168, 1944. a
Liu, J., Bowman, K. W., Schimel, D. S., Parazoo, N. C., Jiang, Z., Lee, M., Bloom, A. A., Wunch, D., Frankenberg, C., Sun, Y., and O'Dell, C. W.: Contrasting carbon cycle responses of the tropical continents to the 2015–2016 El Niño, Science, 358, 6360, https://doi.org/10.1126/science.aam5690, 2017. a
Mackay, D. S., Ewers, B. E., Loranty, M. M., Kruger, E. L., and Samanta, S.: Bayesian analysis of canopy transpiration models: a test of posterior parameter means against measurements, J. Hydrol., 432, 75–83, 2012. a
Mäkelä, A., Berninger, F., and Hari, P.: Optimal control of gas exchange during drought: theoretical analysis, Ann. Bot.London, 77, 461–468, 1996. a
Marquardt, D. W.: An algorithm for leastsquares estimation of nonlinear parameters, J. Soc. Ind. Appl. Math., 11, 431–441, 1963. a, b
McGuire, A. D., Sitch, S., Clein, J. S., Dargaville, R., Esser, G., Foley, J., Heimann, M., Joos, F., Kaplan, J., Kicklighter, D. W., and Meier, R. A.: Carbon balance of the terrestrial biosphere in the twentieth century: Analyses of CO_{2}, climate and land use effects with four processbased ecosystem models, Global Biogeochem. Cy., 15, 183–206, 2001. a
McGuire, A., Wirth, C., Apps, M., Beringer, J., Clein, J., Epstein, H., Kicklighter, D., Bhatti, J., Chapin, F., Groot, B., Efremov, D., Eugster, W., Fukuda, M., Gower, T., Hinzman, L., Huntley, B., Jia, G., Kasischke, E., Melillo, J., Romanovsky, V., Shvidenko, A., Vaganov, E., and Walker, D.: Environmental variation, vegetation distribution, carbon dynamics and water/energy exchange at high latitudes, J. Veg. Sci., 13, 301–314, 2002. a
Medlyn, B. E., Duursma, R. A., Eamus, D., Ellsworth, D. S., Prentice, I. C., Barton, C. V., Crous, K. Y., De Angelis, P., Freeman, M., and Wingate, L.: Reconciling the optimal and empirical approaches to modelling stomatal conductance, Glob. Change Biol., 17, 2134–2144, 2011. a
Monson, R., Turnipseed, A., Sparks, J., Harley, P., ScottDenton, L., Sparks, K., and Huxman, T.: Carbon sequestration in a highelevation, subalpine forest, Glob. Change Biol., 8, 459–478, 2002. a
Monson, R. K., Prater, M. R., Hu, J., Burns, S. P., Sparks, J. P., Sparks, K. L., and ScottDenton, L. E.: Tree species effects on ecosystem wateruse efficiency in a highelevation, subalpine forest, Oecologia, 162, 491–504, 2009. a
Monteith, J. and Unsworth, M.: Principles of environmental physics, Academic Press, USA, 2007. a
Myneni, R. B., Ramakrishna, R., Nemani, R., and Running, S. W.: Estimation of global leaf area index and absorbed PAR using radiative transfer models, IEEE T. Geosci. Remote, 35, 1380–1393, 1997. a
Oleson, K., Lawrence, D., Bonan, G., Drewniak, B., Huang, M., Koven, C., Levis, S., Li, F., Riley, W., Subin, Z., and Swenson, S.: Technical description of version 4.5 of the Community Land Model (CLM), 420 pp., 2013. a
Oleson, K. W., Lawrence, D. M., Gordon, B., Flanner, M. G., Kluzek, E., Peter, J., Levis, S., Swenson, S. C., Thornton, E., Feddema, J., and Heald, C. L.: Technical description of version 4.0 of the Community Land Model (CLM), 2010. a
Pappas, C., Fatichi, S., Leuzinger, S., Wolf, A., and Burlando, P.: Sensitivity analysis of a processbased ecosystem model: Pinpointing parameterization and structural issues, J. Geophys. Res.Biogeo., 118, 505–528, 2013. a
Quaife, T., Lewis, P., De Kauwe, M., Williams, M., Law, B. E., Disney, M., and Bowyer, P.: Assimilating canopy reflectance data into an ecosystem model with an Ensemble Kalman Filter, Remote Sens. Environ., 112, 1347–1364, 2008. a
Reed, D. E., Frank, J. M., Ewers, B. E., and Desai, A. R.: Time dependency of eddy covariance site energy balance, Agr. Forest Meteorol., 249, 467–478, 2018. a
Reichstein, M., Rey, A., Freibauer, A., Tenhunen, J., Valentini, R., Banza, J., Casals, P., Cheng, Y., Grünzweig, J. M., Irvine, J., and Joffre, R.: Modeling temporal and largescale spatial variability of soil respiration from soil water availability, temperature and vegetation productivity indices, Global Biogeochem. Cy., 17, https://doi.org/10.1029/2003GB002035, 2003. a
Ricciuto, D. M., Davis, K. J., and Keller, K.: A Bayesian calibration of a simple carbon cycle model: The role of observations in estimating and reducing uncertainty, Global Biogeochem. Cy., 22, https://doi.org/10.1029/2006GB002908, 2008. a
Rodgers, C. D.: Inverse methods for atmospheric sounding: theory and practice, vol. 2, World scientific, 2000. a, b, c, d
Rogers, A., Medlyn, B. E., Dukes, J. S., Bonan, G., Caemmerer, S., Dietze, M. C., Kattge, J., Leakey, A. D., Mercado, L. M., Niinemets, Ü., and Prentice, I. C.: A roadmap for improving the representation of photosynthesis in Earth system models, New Phytol., 213, 22–42, 2017. a
Running, S. W. and Coughlan, J. C.: A general model of forest ecosystem processes for regional applications I. Hydrologic balance, canopy gas exchange and primary production processes, Ecol. Model., 42, 125–154, 1988. a
Sánchez, J. M., Caselles, V., and Rubio, E. M.: Analysis of the energy balance closure over a FLUXNET boreal forest in Finland, Hydrol. Earth Syst. Sci., 14, 1487–1497, https://doi.org/10.5194/hess1414872010, 2010. a
Schaaf, C. and Wang, Z.: MCD43A4 MODIS/Terra+ Aqua BRDF/Albedo Nadir BRDF Adjusted RefDaily L3 Global 500 m V006, NASA EOSDIS Land Processes DAAC, available at: https://doi.org/10.5067/MODIS/MCD43A4.006, 2015. a
Schimel, D. S.: Terrestrial ecosystems and the carbon cycle, Glob. Change Biol., 1, 77–91, 1995. a
Schulze, E., Kelliher, F. M., Korner, C., Lloyd, J., and Leuning, R.: Relationships among maximum stomatal conductance, ecosystem surface conductance, carbon assimilation rate, and plant nitrogen nutrition: a global ecology scaling exercise, Annu. Rev. Ecol. Syst., 25, 629–662, 1994. a
Simioni, G., Gignoux, J., Le Roux, X., Appé, R., and Benest, D.: Spatial and temporal variations in leaf area index, specific leaf area and leaf nitrogen of two cooccurring savanna tree species, Tree Physiol., 24, 205–216, 2004. a
Sitch, S., Smith, B., Prentice, I. C., Arneth, A., Bondeau, A., Cramer, W., Kaplan, J. O., Levis, S., Lucht, W., Sykes, M. T., and Thonicke, K.: Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model, Glob. Change Biol., 9, 161–185, 2003. a, b
Sitch, S., Huntingford, C., Gedney, N., Levy, P. E., Lomas, M., Piao, S. L., Betts, R., Ciais, P., Cox, P., Friedlingstein, P., and Jones, C. D.: Evaluation of the terrestrial carbon cycle, future plant geography and climatecarbon cycle feedbacks using five Dynamic Global Vegetation Models (DGVMs), Glob. Change Biol., 14, 2015–2039, 2008. a
Sitch, S., Friedlingstein, P., Gruber, N., Jones, S. D., MurrayTortarolo, G., Ahlström, A., Doney, S. C., Graven, H., Heinze, C., Huntingford, C., Levis, S., Levy, P. E., Lomas, M., Poulter, B., Viovy, N., Zaehle, S., Zeng, N., Arneth, A., Bonan, G., Bopp, L., Canadell, J. G., Chevallier, F., Ciais, P., Ellis, R., Gloor, M., Peylin, P., Piao, S. L., Le Quéré, C., Smith, B., Zhu, Z., and Myneni, R.: Recent trends and drivers of regional sources and sinks of carbon dioxide, Biogeosciences, 12, 653–679, https://doi.org/10.5194/bg126532015, 2015. a
Sun, Y., Gu, L., and Dickinson, R. E.: A numerical issue in calculating the coupled carbon and water fluxes in a climate model, J. Geophys. Res., 117, https://doi.org/10.1029/2012JD018059, 2012. a, b, c
Suyker, A. E., Verma, S. B., Burba, G. G., and Arkebauer, T. J.: Gross primary production and ecosystem respiration of irrigated maize and irrigated soybean during a growing season, Agr. Forest Meteorol., 131, 180–190, 2005. a
Tanaka, K., Kosugi, Y., and Nakamura, A.: Impact of leaf physiological characteristics on seasonal variation in CO_{2}, latent and sensible heat exchanges over a tree plantation, Agr. Forest Meteorol., 114, 103–122, 2002. a, b
Tuzet, A., Perrier, A., and Leuning, R.: A coupled model of stomatal conductance, photosynthesis and transpiration, Plant Cell Environ., 26, 1097–1116, 2003. a
van der Tol, C.: SCOPE model, available at: https://github.com/Christiaanvandertol/SCOPE, last access: 28 August 2017.
Van der Tol, C., Dolman, A., Waterloo, M., and Meesters, A.: Optimum vegetation characteristics, assimilation, and transpiration during a dry season: 2. Model evaluation, Water Resour. Res., 44, https://doi.org/10.1029/2007WR006243, 2008a. a
Van der Tol, C., Meesters, A., Dolman, A., and Waterloo, M.: Optimum vegetation characteristics, assimilation, and transpiration during a dry season: 1. Model description, Water Resour. Res., 44, https://doi.org/10.1029/2007WR006241, 2008b. a
van der Tol, C., Verhoef, W., Timmermans, J., Verhoef, A., and Su, Z.: An integrated model of soilcanopy spectral radiances, photosynthesis, fluorescence, temperature and energy balance, Biogeosciences, 6, 3109–3129, https://doi.org/10.5194/bg631092009, 2009. a, b, c, d, e
Verhoef, W., Jia, L., Xiao, Q., and Su, Z.: Unified opticalthermal fourstream radiative transfer theory for homogeneous vegetation canopies, IEEE T. Geosci. Remote, 45, 1808–1822, 2007. a
von Randow, C., Manzi, A. O., Kruijt, B., De Oliveira, P. J., Zanchi, F. B., Silva, R. L., Hodnett, M. G., Gash, J. H. C., Elbers, J. A., Waterloo, M. J., and Cardoso, F. L.: Comparative measurements and seasonal variations in energy and carbon exchange over forest and pasture in South West Amazonia, Theor. Appl. Climatol., 78, 5–26, 2004. a
Wang, Q., Iio, A., Tenhunen, J., and Kakubari, Y.: Annual and seasonal variations in photosynthetic capacity of Fagus crenata along an elevation gradient in the Naeba Mountains, Japan, Tree Physiol., 28, 277–285, 2008. a
Wang, Y.P. and Leuning, R.: A twoleaf model for canopy conductance, photosynthesis and partitioning of available energy I:: Model description and comparison with a multilayered model, Agr. Forest Meteorol., 91, 89–111, 1998. a
Wilson, K., Goldstein, A., Falge, E., Aubinet, M., Baldocchi, D., Berbigier, P., Bernhofer, C., Ceulemans, R., Dolman, H., Field, C., and Grelle, A.: Energy balance closure at FLUXNET sites, Agr. Forest Meteorol., 113, 223–243, 2002. a
Wilson, K. B. and Baldocchi, D. D.: Seasonal and interannual variability of energy fluxes over a broadleaved temperate deciduous forest in North America, Agr. Forest Meteorol., 100, 1–18, 2000. a
Wilson, K. B., Baldocchi, D. D., and Hanson, P. J.: Spatial and seasonal variability of photosynthetic parameters and their relationship to leaf nitrogen in a deciduous forest, Tree Physiol., 20, 565–578, 2000. a, b, c, d, e
Wolf, A., Akshalov, K., Saliendra, N., Johnson, D. A., and Laca, E. A.: Inverse estimation of Vc max, leaf area index, and the BallBerry parameter from carbon and energy fluxes, J. Geophys. Res., 111, 1–18, 2006. a, b, c
Wramneby, A., Smith, B., Zaehle, S., and Sykes, M. T.: Parameter uncertainties in the modelling of vegetation dynamics?effects on tree community structure and ecosystem functioning in European forest biomes, Ecol. Model., 216, 277–290, 2008. a
Wu, X., Luo, Y., Weng, E., White, L., Ma, Y., and Zhou, X.: Conditional inversion to estimate parameters from eddyflux observations, J. Plant Ecol., 2, 55–68, 2009. a
Wullschleger, S. D.: Biochemical limitations to carbon assimilation in C3 plants – a retrospective analysis of the A/Ci curves from 109 species, J. Exp. Bot., 44, 907–920, 1993. a, b
Xu, L. and Baldocchi, D. D.: Seasonal trends in photosynthetic parameters and stomatal conductance of blue oak (Quercus douglasii) under prolonged summer drought and high temperature, Tree Physiol., 23, 865–877, 2003. a, b
Xu, T., White, L., Hui, D., and Luo, Y.: Probabilistic inversion of a terrestrial ecosystem model: Analysis of uncertainty in parameter estimation and model prediction, Global Biogeochem. Cy., 20, https://doi.org/10.1029/2005GB002468, 2006. a
Zaehle, S., Sitch, S., Smith, B., and Hatterman, F.: Effects of parameter uncertainties on the modeling of terrestrial biosphere dynamics, Global Biogeochem. Cy., 19, https://doi.org/10.1029/2004GB002395, 2005. a
 Abstract
 Introduction
 SCOPE model
 Formulation of inverse problem
 Linearization of the forward model
 Iterative retrieval algorithm setup
 Results for implementing the inversion framework in SCOPE
 Discussion and conclusion
 Code and data availability
 Appendix A: Modeling photosynthesis in SCOPE
 Appendix B: Temperature dependence of photosynthetic parameters
 Appendix C: Derivation of iterative retrieval algorithm
 Author contributions
 Competing interests
 Acknowledgements
 References
 Supplement
 Abstract
 Introduction
 SCOPE model
 Formulation of inverse problem
 Linearization of the forward model
 Iterative retrieval algorithm setup
 Results for implementing the inversion framework in SCOPE
 Discussion and conclusion
 Code and data availability
 Appendix A: Modeling photosynthesis in SCOPE
 Appendix B: Temperature dependence of photosynthetic parameters
 Appendix C: Derivation of iterative retrieval algorithm
 Author contributions
 Competing interests
 Acknowledgements
 References
 Supplement