Optimal inverse estimation of ecosystem parameters from observations of carbon and energy fluxes

Canopy structural and leaf photosynthesis parameterizations such as maximum carboxylation capacity (Vcmax), slope of the Ball–Berry stomatal conductance model (BBslope) 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 Vcmax, BBslope 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 C3 and C4 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 (0.95>R2 > 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.

Abstract. 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 (0.95 > R 2 > 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.

Introduction
Terrestrial ecosystems play a very important role in regulating the carbon exchange over land surfaces (Schimel, 1995;Falkowski et al., 2000). 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(Baldocchi et al., , 2001Sitch et al., 2003Sitch et al., , 2008. Terrestrial ecosystem models have been used to study the carbon and water fluxes (McGuire et al., 2001;Sitch et al., 2003;Cramer et al., 2001;Kucharik et al., 2000); however, there are large uncertainties in fluxes associated with poorly quantified model parameters (Wramneby et al., 2008;Pappas et al., 2013;Knorr and Heimann, 2001;Zaehle et al., 2005;Rogers et al., 2017). Some of these parameters have temporal and spatial variability and are hard to measure directly over large scales (Simioni et al., 2004;Dutta et al., 2017). 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 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 (van der Tol et al., 2009;Running and Coughlan, 1988;Oleson et al., 2010). Detailed canopy models are often resolved into multiple layers, thus providing a better treatment of radiation regime and energy balance across the canopy (van der Tol et al., 2009;Wang and Leuning, 1998;Dai et al., 2004). At the heart of these models lies the leaf-level 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 process-based 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 enzyme-catalyzed reaction of the substrates CO 2 or O 2 with ribulose-1,5-bisphosphate, representing the enzymelimited photosynthesis rate (Farquhar et al., 1980). Other rate-limiting 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;Wullschleger, 1993;Beerling and Quick, 1995;Tanaka et al., 2002). 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 leaf-level gas exchange measurements and generated A − C i curves (Wullschleger, 1993;Tanaka et al., 2002;Xu and Baldocchi, 2003). LAI can be estimated from destructive and nondestructive optical methods (Myneni et al., 1997;Dutta et al., 2017;Chen 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 process-based models using observations of carbon and energy fluxes could thus yield these key ecosystem parameters. Process-based 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 solar-induced 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 least-squares 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 radiation-driven photosynthesis, respiration and energy balance using a two-component (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;Xu et al., 2006;Ricciuto et al., 2008;Mackay et al., 2012) 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 half-yearly 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 sun-induced chlorophyll fluorescence (SIF). SIF is chlorophyll re-emission during photosynthesis, acts as a direct probe into photosynthesis measurable from space and is strongly correlated with flux-based GPP estimates at canopy to ecosystem scales (Frankenberg et al., 2011;Flexas et al., 2002). Thus, the SCOPE-based inversion approach has the flexibility and advantage of incorporating tower-based 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: 1. implementation of photosynthesis model and its temperature dependencies consistent with a well-accepted major Earth system model (Community Land Model CLM4.5) in SCOPE; 2. development of a Bayesian nonlinear inversion framework using SCOPE to estimate ecosystem parameters using eddy covariance flux observations; 3. 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.

SCOPE model
The SCOPE model (van der Tol et al., 2009) is an integrated 1-D 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 top-of-canopy radiance spectrum of fluorescence from leaf-level chlorophyll fluorescence. SCOPE resolves top-of-canopy 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.

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) Oleson et al., 2013) include updates to the canopy radiation scheme and canopy scaling of leaf processes, colimitations on photosynthesis, revisions to photosynthetic parameters , 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 multi-layer canopy modeling framework with coupled photosynthesis (Farquhar et al., 1980) and Ball-Berry stomatal conductance models similar to the SCOPE framework. Figure 1. Temperature response functions of V cmax for C 3 (a) and C 4 (b) plants. For the C 3 species the mean and ±1σ variability (shown as broken black lines) in the net temperature response are computed using data presented in Leuning (2002), shown in the left panel. The temperature range corresponding to maximum V cmax response for both the C 3 and C 4 pathways is between 30 and 40 • C. The overall temperature response from the previous version (V1.70) is shown as dashed brown line.
The main inconsistencies of SCOPE (V1.70) with the CLM4.5 parameterizations are as follows: 1. 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.
2. 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.
3. The net assimilation, internal CO 2 concentration and stomatal conductance (A-C 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: 1. 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 ).
2. The light limited photosynthesis rate for C 4 is given by Eq. (A3).
3. The temperature dependence of photosynthetic parameters  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 (Leuning, 2002;Bernacchi et al., 2001) and C 4 photosynthesis is represented by Eqs. (B1)-(B5).
4. Finally we also incorporate a new simplified implementation of A-C 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. 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 photosynthesis-pathway-dependent 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 . A number of analyses were performed to study the differences in net response of canopy-level 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 GPP = GPP new GPP 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 2-D histograms, of which only the mean (f GPP mean ) is represented in Fig. 2.

Comparison of current and previous photosynthesis implementations in SCOPE
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.

Formulation of inverse problem
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 F() : X −→ 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 F(), 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 (X, p ⊂ X and X = X ∪ 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 vectorX can be obtained by a generalized inverse method R represented as follows: wherep 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 catch-all for anything else that is used in the retrieval method, which also includes the convergence criteria. ) simulations as a function of PAR, canopy temperature and VPD for the C 3 species. The Missouri Ozark fluxnet site comprising of deciduous broadleaf forests for the year 2009 is used and the SCOPE simulations are driven with identical forcings and parameters for both the new and old simulations, the only difference being the implementation of photosynthesis and its temperature dependence (see text). The left column shows f GPP as a function of PAR and canopy temperature with data points in the low VPD range of 10-15 hPa (a) and high VPD range of 25-30 hPa (b). The right column shows f GPP as a function of PAR and VPD with data points in the low temperature range of 10-15 • C (c) and moderate temperature range of 25-30 • C (d).

Linearization of the forward model
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 first-order Taylor expansion for the forward model where X l is an arbitrary linearization point, and δF δX is the partial derivative or Jacobian at the point X = X l .

Iterative retrieval algorithm setup
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 vectorX (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: BB slope = 7 and LAI = 4. Second, third and fourth rows from the top shows the diurnal variability in the gradient 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 which constitutes the Jacobian matrix for the inversions. It can be observed that the Jacobian matrix is nonlinear with maximum values near the midday period. Our retrieval framework uses concatenated 3-day GPP and LE fluxes (modeled and observed) and their gradients successively within a 3-day window.
where K i is the Jacobi matrix. A brief derivation of Eq. (6) is presented in Appendix C; for a more in-depth treatment the reader is referred to Rodgers (2000).

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 least-squares 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 S −1 a (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:

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 Mead-1 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 perturba- tions δ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 leaf-level 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 half-hourly 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 long-term (seasonal) retrieval of im-portant 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 N-day contiguous moving windows over the entire growing season (Fig. 4). We extend the 1-day 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 N-day 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 3-day 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.

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 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 vectorX can be computed as follows: The reduction in error is defined as follows: where S jj and S 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 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.

Results for implementing the inversion framework in SCOPE
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.

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.
1. In constructing the observation vector Y we apply a time-of-the-day filter (e.g. data between 09:00 and 16:00 LT and so on) for the initial forward SCOPE model.
2. For computing the Jacobians, a PAR-based 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).
3. 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.

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 Mead-1 site

Site description
The Nebraska Mead-1 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.

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 off-diagonal elements. Figure 7 shows the retrievals of parameters V cmax , BB slope and LAI. 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 ζ V cmax and ζ LAI values and moderate to no seasonality for the ζ BB 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 off-diagonal elements). The error correlation is given by ρ x,y = COV(x,y) σ x σ y and should be considered when interpreting covariations of retrieved parameters as the nature and mag-  nitude 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 ρ BB slope,LAI and ρ V cmax,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 ρ V cmax ,BB 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 one-to-one 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 ( 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 half-hourly 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 late-summer 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).

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 Mead-1, 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 ( = |Observed − 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 ζ V 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 Corn-C 4 site using flux observations only shows it is similar for ρ V cmax,LAI and ρ BB slope,LAI (negative and positive respectively) and different for ρ V cmax ,BB slope .
The correlations ρ V cmax ,BB slope and ρ V cmax,LAI are both negative, indicating counteracting effects of these variables in the retrieval. In comparison, ρ BB slope,LAI is positive for the retrievals, indicating in-sync behavior. The error correlations ρ V cmax,LAI are highest among the three in the middle of the growing season.  and χ 2 LE−unopt = 32 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.

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 late-summer drought. The 92 D. Dutta et al.: Optimal inverse estimation of ecosystem parameters 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 post-drought 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 2 GPP−opt = 0.7, R 2 GPP−unopt = 0.2, R 2 LE−opt = 0.5 and R 2 LE−unopt = 0.1. For the year 2007 it is also found that the posterior error reduction ζ V 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 ρ V cmax ,BB slope and ρ V cmax,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 in-season drought dynamics, which are not possible using traditional one-step 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. 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 half-hourly 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.  Panel (a) shows the ζ V cmax , ζ BB slope and ζ LAI for the entire growing season and (b) shows the correlation coefficients (normalized offdiagonal elements of posterior error covariance matrix) among these variables. Both the ζ and correlation coefficients are computed using the final Jacobian matrix at the end of each retrieval window.

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 fluctu-ations. 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 ρ V cmax ,BB slope and ρ V cmax,LAI are both negative and ρ BB slope,LAI is positive. There is a sharp discontinuity around DOY 250-260 in terms of ζ BB slope and ρ V cmax,LAI , ρ V cmax ,BB 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 2 GPP−opt = 0.85, R 2 GPP−unopt = 0.79, R 2 LE−opt = 0.52 and R 2 LE−unopt = 0.49. The corresponding error measures are χ 2 GPP−opt = 38 971, χ 2 GPP−unopt = 22 085, χ 2 LE−opt = 63 865 and χ 2 LE−unopt = 68 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 near-constant 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.

Discussion and conclusion
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;, 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  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 . These changes may be mostly attributed to substantial in-season changes in the fraction of total N allocated to rubisco as well as changes in leaf mass per area . 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 sub-alpine 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 % (Wilson et al., 2002;von Randow et al., 2004;Sánchez et al., 2010) 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 (Van der Tol et al., 2008b, a;Mäkelä et al., 1996) 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 (Medlyn et al., 2011;Katul et al., 2009).
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 space-based (Frankenberg et al., 2011Guanter et al., 2014) and ground-based 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.
Code and data availability. The authors thank the AmeriFlux team for making the eddy-covariance 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/fluxnet2015-dataset/, 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: for C 4 species.
The light-limited rate of carboxylation (governed by the capacity to regenerate RuBP) A j is given by the following: 4(C i +2 * ) , for C 3 species, α(4.6φ), for C 4 species.
Finally the product-limited carboxylation rate for C 3 plants and the PEP-carboxylase-limited rate of carboxylation for the C 4 plants A p is given by the following: A p = 3T p , for C 3 species, k p C i , for C 4 species.
For the above Eqs. A2, A3 and A4, we have the assimilation rates A c,j,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.
I PSII = 0.5 PSII (4.6φ) The term PSII in Eq. (A6) is the quantum yield of photosystem II and 0.5 represents half-electron transfer to each of the photosystems I and II. The overall gross photosynthesis rate is computed as a colimitation (Collatz et al., 1991a;Collatz et al., 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 light-limited and enzymeor product-limiting 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 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.
B1 A-g 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.
C i = C s − 1.6 A g s (from Fick's law), (B6) g s = g 0 + BB slope Ar h C s (Ball-Berry model).
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 min Ci 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 (Sun et al., 2012;Ivanov et al., 2008). The updated C i is further used in the A-g s -C i until convergence.

98
D. Dutta et al.: Optimal inverse estimation of ecosystem parameters Table B1. C 3 temperature dependence functional parameters.
Parameter Q 10 (-) s 1 (K −1 ) s 2 (K) s 3 (K −1 ) s 4 (K) Author contributions. 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.
Competing interests. The authors declare that they have no conflict of interest.