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.

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 CO2
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 inter-annual 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 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 CO2 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 (Vcmax), which is directly related to the
concentration of the enzyme rubisco. Vcmax is a key parameter
in the Michaelis–Menten kinetics for an enzyme-catalyzed reaction of the
substrates CO2 or O2 with ribulose-1,5-bisphosphate,
representing the enzyme-limited photosynthesis rate
(Farquhar et al., 1980). Other rate-limiting photosynthesis
parameters such as maximum electron transport rate (Jmax) are
generally parameterized with respect to Vcmax. The Ball–Berry
equation calculates the stomatal conductance (gs) for water vapor
as a function of net assimilation, relative humidity, leaf surface
CO2 concentration, minimum conductance and a proportionality
constant called the Ball–Berry slope (BBslope)
(Ball et al., 1987; Beerling and Quick, 1995; Tanaka et al., 2002; Wullschleger, 1993).
The BBslope 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 Vcmax and
BBslope can be determined experimentally from leaf-level gas
exchange measurements and generated A−Ci 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 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 (Vcmax,
BBslope 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
(Vcmax, BBslope, 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 Vcmax and respiration rate and energy difference
(H − LE) to fit BBslope. 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
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
(Flexas et al., 2002; Frankenberg et al., 2011). 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:

implementation of photosynthesis model and its temperature dependencies
consistent with a well-accepted 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 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.

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 community-developed 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
multi-layer canopy modeling framework with coupled photosynthesis
(Farquhar et al., 1980) and Ball–Berry stomatal conductance models
similar to the SCOPE framework.

Figure 1Temperature response functions of Vcmax for
C3(a) and C4(b) plants. For the C3 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 Vcmax response for both the
C3 and C4 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:

Similar, generic temperature response functions are implemented for
both C3 and C4 species, excepting Vcmax, and further it
uses a Q10-based exponential function with the same functional parameters
for computing the temperature response of the various photosynthetic parameters.

There is no Jmax (maximum potential electron transport
rate, ETR) or its temperature dependence in the computation of light limited
C3 photosynthesis rate.

The net assimilation, internal CO2 concentration and stomatal
conductance (A-Ci-gs) 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 Vcmax 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
Vcmax. Overall, the major new updates made to the model
(biochemical module) are as follows:

The electron limited photosynthesis rate Aj 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
(IPSII) and the maximum potential ETR (Jmax).

Finally we also incorporate a new simplified implementation
of A-Ci-gs 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 Vcmax for both C3 (left) and C4
(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 C3 and the Q10 functions for the C4 pathways.
We have also shown the temperature dependence of Vcmax 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 C3 and C4 species. It can be observed that
the difference in temperature response is more for C4, 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).

Figure 2Figure showing the ratio
of old and new SCOPE GPP (fGPP=GPPnewGPPold)
simulations as a function of PAR, canopy temperature and VPD for the C3
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 fGPP 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 fGPP 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).

A number of analyses were performed to study the differences in net response
of canopy-level carbon and energy fluxes for both C3 and C4 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 fGPP=GPPnewGPPold.
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 C3 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 fGPP 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 (fGPPmean) is represented in Fig. 2.

We find that over the larger parts of the domain of random variables,
fGPP 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
C3 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 Vcmax
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), fGPP 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 Vcmax 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 Vcmax response in that temperature range as indicated
in Fig. 1. The results for a similar comparative study with C4 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−Ci 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:

(1)Y=F(X′)+ε,

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 ℱ(), which is complex and moderately
nonlinear, representing a range of physics and canopy physiological
processes. We can further represent our forward problem as follows:

(2)Y=F(X;p)+ε,

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
vector X^ can be obtained by a generalized inverse method
ℛ represented as follows:

(3)X^=R(Y,p^,Xa,c),

where p^ represents the best estimate of the forward function
parameters. The parameters Xa and c represent the
parameters that do not appear in the forward function but they do affect the
retrieval and are associated with uncertainties. Xa
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.

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

(4)F(X;p)=F(X;p)X=Xl+δFδXX=Xl(X-Xl),

where Xl is an arbitrary linearization point, and
δFδX is the partial derivative or Jacobian at
the point X=Xl.

Figure 3Diurnal variability of GPP, LE, H and SIF from SCOPE model
simulations (top row) for a typical day in the growing season (3 August 2010)
for C4 corn using data from the Nebraska Mead-1 flux tower site with
parameter values Vcmax=50µmols m−2 s−1,
BBslope=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 δVcmax=5µmols m−2 s−1,
δBBslope=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.

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 (Xa) with a Gaussian uncertainty in the prior state
(Sa) (Rodgers, 2000) can be represented as
follows:

-2lnP(X|Y)=[Y-F(X)]TSε-1[Y-F(X)](5)+[X-Xa]TSa-1[X-Xa]+c′,

where c′ is a constant. Our aim is to find the best estimate of the state
vector 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:

Xi+1=Xi+(Sa-1+KiTSε-1Ki)-1[KiTSε-1[Y-F(Xi)](6)-Sa-1[Xi-Xa]],

where Ki 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).

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 least-squares problem:

(7)Xi+1=Xi+(KKT+γiD)-1KT[Y-F(Xi)],

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
Sa-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:

Xi+1=Xi+[(1+γ)Sa-1+KiTSε-1Ki]-1(8){KiTSε-1[Y-F(Xi)]-Sa-1[Xi-Xa]}.

Figure 4Illustration of moving window inversion retrieval setup. The bottom
left part illustrates the annual ecosystem time series flux variables used
for driving the SCOPE model. A 3-day time window is selected for each
retrieval in the yearly growing season and a time filter is implemented for
concatenating the measurement vector (flux tower observations, with each
color representing a different observation variable) of length m in the
retrieval windows. The top right shows the vector and matrix setup and the
linearization of the forward model. ΔY represents the
difference between concatenated observation and modeled vector and ΔX represents the corresponding change in the state vector comprising
of n variables (parameters). The bottom right shows the retrieved model
parameters after implementing the moving window approach.

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 C4 corn
using data from the Nebraska Mead-1 flux tower site with parameter values
Vcmax=50µmols m−2 s−1,
BBslope=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
δVcmax=5µmols m−2 s−1,
δBBslope=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, BBslope 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 Vcmax 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 Ym×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 (Xa[n×1]) as well as the prior error covariance matrix
(Sa[n×n]). As mentioned previously, the
Jacobian matrix Km×n is computed numerically by a small
perturbation to the value of the state vector Xi+δ (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 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 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 N-day 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 (Vcmax, BBslope and LAI)
towards the change in environmental drivers.

Figure 5Figure
showing the seasonal variability of red, NIR MODIS nadir reflectance and NDVI
for the Nebraska Mead-1 and Missouri Ozark Site. This dataset is used in SCOPE
model inversions in conjunction with fluxes of carbon and water for retrieval
of ecosystem parameters.

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
[Y-F(Xi)]TSε-1[Y-F(Xi)]+[X-Xa]TSa-1[X-Xa].
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 X^ can be computed as follows:

(9)S=[Sa-1+KiTSε-1Ki]-1.

The reduction in error is defined as follows:

(10)ζi=1-SjjSajj0.5,

where Sjj and Sajj 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.

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
Vcmax, BBslope 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 time-of-the-day
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 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).

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.

Figure 6Figure showing the diurnal and seasonal variability of important
environmental and meteorological forcings together with the tower observed
fluxes of carbon and energy used in SCOPE model inversions for the Nebraska
Mead flux tower site. The variables in (a) and (b) are used
as inputs to the SCOPE model and the variables in (c) are used as a
target in a moving window retrieval approach.

6.3 Retrieval results for the Nebraska Mead-1 site

6.3.1 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 (C4, 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 Vcmax, BBslope
and LAI during this entire growing season.

Figure 7Figure showing the seasonal variability in retrieved parameter
values of Vcmax, BBslope and LAI for the Nebraska
Mead-1 site using a 3-day moving window inversion approach for the year 2010.
The actual points in the time series (grey lines) of the GPP and LE fluxes
used as the target observations (Y) for the moving window inversion
approach are shown in the background. The figure shows the comparison of the
retrieved parameters using only GPP and LE fluxes (shown as dashed lines) as
well as using a combination of fluxes and MODIS reflectance (shown as solid
lines). The results show reasonable trends in the retrieved parameters along
with their sensitivity to GPP and LE fluxes across the 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 off-diagonal
elements. Figure 7 shows the retrievals of parameters
Vcmax, BBslope 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 (Δ=|Observed-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
Vcmax shows a similar seasonality to GPP. There is also some
variability in retrieved BBslope, 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 Vcmax.

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
Vcmax. 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
Vcmax, BBslope 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 BBslope
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 BBslope 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 BBslope, even though it would not
represent the gas exchange through the stomata. From the inversion results it
is clear that all three parameters Vcmax, BBslope
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 8Figure showing the
seasonal variability of the posterior error reduction (ζ) and
correlation coefficient of the retrieved parameter values of
Vcmax, BBslope and LAI for the Nebraska Mead-1
site using a 3-day moving window inversion approach for the year 2010. Panel (a) shows the ζVcmax,
ζBBslope and ζLAI for the entire
growing season and (b) shows the correlation coefficients
(normalized off-diagonal elements of posterior error covariance matrix) among
these variables. The results of retrievals using only GPP and LE fluxes are
shown as dashed lines and the results using a combination of fluxes and MODIS
reflectance are shown as solid lines. Both the ζ and correlation
coefficients are computed using the final Jacobian matrix at the end of each
retrieval window.

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 ζVcmax and
ζLAI values and moderate to no seasonality for the
ζBBslope. 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
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 ρBBslope,LAI and
ρVcmax,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
ρVcmax,BBslope 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 (R2=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 (χ2=∑i=0k((yi-F(xi))/σi)2) for the
prior (unoptimized) and posterior (optimized) simulations. The corresponding
values are χGPP-opt2=9382, χGPP-unopt2=10728,
χLE-opt2=19235 and χLE-unopt2=20554.

Figure 9Figure showing the
improvement in diurnal and seasonal variability in modeling the GPP and LE
fluxes with optimized parameters over prior values using SCOPE for the
Nebraska Mead-1 site for the year 2010. The figure also shows the one-to-one
comparison (indicated by the black line) with the observed flux tower values. The
optimization of the photosynthetic parameters improves the accuracy of
computing the carbon and water fluxes as indicated by the R2 value and the
equation of the regression line.

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 C3 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 Vcmax fits
reflect the stress factor as well, which is usually applied to downscale
the physiological Vcmax during environmental stress. We focus
on the retrieval of the parameters Vcmax, BBslope
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).

Figure 10Figure showing the seasonal variability in retrieved parameter
values of Vcmax, BBslope and LAI for the Missouri
Ozark site using a 3-day moving window inversion approach for the year 2009.
The actual points in the time series (grey lines) of the GPP and LE fluxes
used as the target observations (Y) for the moving window inversion
approach are shown in the background. The figure shows the comparison of the
retrieved parameters using only GPP and LE fluxes (shown as dashed lines) as
well as using a combination of fluxes and MODIS reflectance (shown as solid
lines). The results show reasonable trends in the retrieved parameters along
with their sensitivity to GPP and LE fluxes across the growing season.

Figure 11Figure showing the seasonal variability of the posterior error
reduction (ζ) and correlation coefficient of the retrieved parameter
values of Vcmax, BBslope and LAI for the Missouri
Ozark site using a 3-day moving window inversion approach for the year
2009. Panel (a) shows the ζVcmax,
ζBBslope and ζLAI for the
entire growing season and (b) shows the correlation coefficients
(normalized off-diagonal elements of posterior error covariance matrix) among
these variables. The results of retrievals using only GPP and LE fluxes are
shown as dashed lines and the results using a combination of fluxes and MODIS
reflectance are shown as solid lines. Both the ζ and correlation
coefficients are computed using the final Jacobian matrix at the end of each
retrieval window.

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 Mead-1,
Fig. 7) shows the comparison of parameter retrievals
Vcmax, BBslope 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
Vcmax, BBslope and LAI, following the patterns in
GPP and LE. The beginning (and end) of the growing season shows increasing
(decreasing) trends in Vcmax 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 Vcmax, which is specifically
evident around DOYs 140–150 and 200–250. The sharp increase in retrieved
Vcmax 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
(BBslope); this may again correspond to the large precipitation
(see Fig. 4) events (e.g. overestimation of BBslope) 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 12Figure showing the seasonal variability in retrieved parameter
values of Vcmax, BBslope and LAI for the Missouri
Ozark site using a 3-day moving window inversion approach for the year 2007.
The actual points in the time series (grey lines) of the GPP and LE fluxes
used as the target observations (Y) for the moving window inversion
approach are shown in the background. The figure shows the comparison of the
retrieved parameters using only GPP and LE fluxes (shown as dashed lines) as
well as using a combination of fluxes and MODIS reflectance (shown as solid
lines). The results show reasonable trends in the retrieved parameters along
with their sensitivity to GPP and LE fluxes across the growing season.

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 BBslope. We
find that ζVcmax and ζLAI have the same
trend and seasonality as that of retrieved Vcmax and LAI
respectively. The evolution of the retrieval error correlations is again
interesting and comparison with the Mead Corn-C4 site using flux
observations only shows it is similar for ρVcmax,LAI
and ρBBslope,LAI (negative and positive respectively) and
different for ρVcmax,BBslope.

The correlations ρVcmax,BBslope and
ρVcmax,LAI are both negative, indicating counteracting
effects of these variables in the retrieval. In comparison,
ρBBslope,LAI is positive for the retrievals, indicating
in-sync behavior. The error correlations ρVcmax,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) (R2=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 χGPP-opt2=12975,
χGPP-unopt2=30070, χLE-opt2=17260 and
χLE-unopt2=32283 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
Vcmax, BBslope 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 midseason drought around DOY 230–250 is well captured by decreases in photosynthesis (Vcmax)
and stomatal conductance (BBslope). These parameters again
increase from around DOY 260 after the drought recovery and also correlate
well with the flux observations. However, the Vcmax 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 Vcmax
with a stress factor between 0 and 1. Our retrieved effective Vcmax
is the product of a stress factor and the physiological Vcmax,
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 Vcmax) much
better than just the flux observations, which is clearly evident from
the Vcmax retrievals around DOY 200. This large change in
Vcmax (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 Vcmax 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
RGPP-opt2=0.7, RGPP-unopt2=0.2,
RLE-opt2=0.5 and RLE-unopt2=0.1. For the year
2007 it is also found that the posterior error reduction
ζVcmax 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 ρVcmax,BBslope and ρVcmax,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 down-regulation 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
Vcmax variations due to changes in the average canopy
temperatures well. There is a clear difference between the retrieved
Vcmax 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.

Figure 13Figure showing the seasonal variability in retrieved parameter
values of Vcmax, BBslope and LAI for the Niwot
Ridge site using a 3-day moving window inversion approach for the year 2010.
The actual points in the time series (grey lines) of the GPP and LE fluxes
used as the target observations (Y) for the moving window inversion
approach are shown in the background. The results show reasonable trends in
the retrieved parameters along with their sensitivity to GPP and LE fluxes
across the growing season.

Figure 14Figure
showing the seasonal variability of the posterior error reduction (ζ)
and correlation coefficient of the retrieved parameter values of
Vcmax, BBslope and LAI for the Niwot Ridge site
using a 3-day moving window inversion approach for the year 2010. Panel (a) shows the ζVcmax, ζBBslope
and ζLAI for the entire growing season and (b) shows the correlation coefficients (normalized off-diagonal 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.

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 C3 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
Vcmax and BBslope 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
Vcmax, BBslope 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 Vcmax 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
BBslope 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 BBslope,
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 BBslope and Vcmax 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 C3
site. We find ρVcmax,BBslope and
ρVcmax,LAI are both negative and
ρBBslope,LAI is positive. There is a sharp discontinuity
around DOY 250–260 in terms of ζBBslope and
ρVcmax,LAI, ρVcmax,BBslope, 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 RGPP-opt2=0.85,
RGPP-unopt2=0.79, RLE-opt2=0.52 and
RLE-unopt2=0.49. The corresponding error measures are
χGPP-opt2=38971, χGPP-unopt2=22085,
χLE-opt2=63865 and χLE-unopt2=68463.
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 Vcmax). 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
Vcmax). 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 Vcmax, 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 Vcmax and
BBslope 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 Vcmax and BBslope (to some
extent). Previously, many studies have reported measured Vcmax
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
Vcmax 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 in-season 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
Vcmax by a stress factor [0, 1]. The effect of reductions in Vcmax are a reduction
in assimilation, which also suppresses transpiration. Thus, we are fitting an
effectiveVcmax 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 BBslope (15 to 25) values
(Wolf et al., 2006). The SCOPE model handles both the C3 and C4
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 % (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
BBslope(Mäkelä et al., 1996; Van der Tol et al., 2008a, b) or
Vcmax(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 Vcmax). Global time series of SIF
observations, which provide a direct probe into photosynthetic machinery, are
now available from space-based
(Frankenberg et al., 2011, 2014; Guanter 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 Vcmax. 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 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 C3 and C4 plant
species and is represented as follows:

The net photosynthesis (assimilation) after accounting for respiration (Rd) is given as follows:

(A1)An=min(Ac,Aj,Ap)-Rd.

Further, the rate limiting steps are represented as follows.
The RuBP carboxylase (rubisco) limited rate of carboxylation Ac is given by they following:

(A2)Ac=Vcmax(Ci-Γ*)Ci+Kc(1+Oi/Ko),forC3species,Vcmax,forC4species.

The light-limited rate of carboxylation (governed by the capacity to regenerate RuBP) Aj is given by the following:

(A3)Aj=J(Ci-Γ*)4(Ci+2Γ*),forC3species,α(4.6ϕ),forC4species.

Finally the product-limited carboxylation rate for C3 plants and the
PEP-carboxylase-limited rate of carboxylation for the C4 plants Ap is
given by the following:

(A4)Ap=3Tp,forC3species,kpCi,forC4species.

For the above Eqs. A2, A3 and A4, we have
the assimilation rates Ac,j,p (in the units of
µmols m−2 s−1), Ci is the internal
CO2 concentration of the leaf (units of ppm) and
Vcmax is the maximum rate of carboxylation. For the C3
species, Kc and Ko are the Michaelis–Menten
constants for CO2 and O2 respectively (units of
µmols m−2 s−1), Γ* is the CO2
compensation point (units are ppm), J is the potential electron transport
rate (units of µmols m−2 s−1) and Tp is the
triose phosphate utilization rate. For the C4 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
CO2 mol−1 photon), and kp is the initial slope of
the C4CO2 response curve.

For the C3 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:

(A5)ΘPSIIJ2-(IPSII+Jmax)J+IPSIIJmax=0,

where Jmax is the maximum electron transport rate (µmols m−2 s−1),
IPSII is the light used in photosystem II (µmols m−2 s−1)
which is given by Eq. (A6) and ΘPSII is a
curvature parameter.

(A6)IPSII=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, 1992) and is computed as the
smaller root of the following equations:

ΘcjAi2-(Ac+Aj)Ai+AcAj=0,(A7)ΘipA2-(Ai+Ap)A+AiAp=0.

The parameters Θcj and Θip control the
smoothness of the light response curve between light-limited and enzyme- or
product-limiting rates. The values of the different parameters
at optimum temperature (mostly as a function of Vcmax25, 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 C3
and C4 pathways described in the previous section and shown in Table A1 have temperature-dependent 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 C3 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:

f(Tv)=expHaRT0(1-T0Tv),(B1)fH(Tv)=1+expSvT0-HdRT01+expSvTv-HdRTv,

where Ha is the activation energy, Hd is the
deactivation energy, Sv is the entropy term, T0 is the optimum
temperature and Tv is the leaf temperature. The functional relationship of
the different photosynthetic parameters in the C3 pathway are as follows:

The Q10 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 C3 and C4 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.

In Eq. (B7), BBslope represents the Ball–Berry
slope, rh the relative humidity and g0 the Ball–Berry intercept. In the
absence of an initial specification of Ci, we make the assumption
that g0=0 in Eq. (B7); then, combining Eqs. (B6), (B7), the initial estimate of Ci is
given as follows:

(B8)Ci=maxfCiminCs,Cs-1.6CsrhBBslope,

where fCimin is the assumed minimum fractional leaf
boundary layer CO2 (assumed to be 0.3 for C3 and 0.1 for C4
species). This initial estimate of Ci 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 CO2 concentration
(Ivanov et al., 2008; Sun et al., 2012). The updated Ci is
further used in the A-gs-Ci until convergence.

For deriving the maximum probability state X (X^) we
equate the derivative of the Eq. (5) to zero to obtain the following:

∇X{-2lnP(X|Y)}=-[∇XF(X)]TSε-1[Y-F(X)](C1)+Sa-1[X-Xa]=0.

It can be noted here that the gradient ∇X of the above
vector-valued function is a matrix-valued function and the Jacobian matrix is
represented as K(X)=∇XF(X), which
results in the following implicit equation for X^:

(C2)-KT(X)Sε-1[Y-F(X)]+Sa-1[X-Xa]=0.

We have to now use any general root-finding 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:

(C3)Xi+1=Xi-[∇XG(Xi)]-1G(Xi).

For our problem we can assume the derivative of the cost-function
G(X) to be the left-hand side of Eq. (C1); therefore the
gradient of G(X) (∇G) also known as the Hessian is given by the following:

∇XG(X)=(C4)Sa-1+KTSε-1K-[∇XKT]Sε-1[Y-F(X)].

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 post-multiplied by the factor
Sε-1[Y-F(X)]. The third
term in the right-hand 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 small-residual
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):

Xi+1=Xi+(Sa-1+KiTSε-1Ki)-1[KiTSε-1[Y-F(Xi)](C5)-Sa-1[Xi-Xa]],

where Ki=K(Xi). We can substitute
F(X) from Eq. (4) in Eq. (C2) to
get the following:

-KT(X^)Sε-1[Y-F(X)X=Xl+∇XF(X)X=Xl(X-Xl)](C6)+Sa-1[X^-Xa]=0.

Again, representing Kl=∇XF(X)X=Xl,
Fl=F(X)X=Xl we can further
simplify and rearrange Eq. (C6) as follows:

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.

AmeriFlux site Missouri Ozark (US-MOz) 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 Science-Science Focus Area; ORNL is managed by UT-Battelle, LLC,
for DOE under contract DE-AC05-00OR22725.

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.

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

Canopy structural and leaf photosynthesis parameterizations are often fixed over time in Earth system models and represent large sources of uncertainty in predictions of carbon and water fluxes. We develop a moving window nonlinear optimal parameter inversion framework using constraining flux and satellite reflectance observations. The results demonstrate the applicability of the approach for error reduction and capturing the seasonal variability of key ecosystem parameters.

Canopy structural and leaf photosynthesis parameterizations are often fixed over time in Earth...