2.1 BETHY-SCOPE

BETHY-SCOPE is a coupling of the existing models BETHY (Biosphere Energy Transfer Hydrology) (Knorr, 2000) and SCOPE (Soil Canopy Observation, Photosynthesis and Energy fluxes; van der Tol et al., 2009, 2014) and builds upon the developments by Koffi et al. (2015) and Norton et al. (2018). Here, some model developments have been made that distinguish this version of BETHY-SCOPE from that used in Norton et al. (2018); hence we refer to this version as BETHY-SCOPE v1.1. The coupling of BETHY and SCOPE enables spatially explicit global simulations of GPP and SIF that are dependent on plant functional type (PFT).

BETHY is a process-based terrestrial biosphere model, which is a key element of the CCDAS (Rayner et al., 2005; Scholze et al., 2007). Full model description details can be found elsewhere (e.g., Rayner et al., 2005; Scholze et al., 2007; Knorr et al., 2010). Briefly, BETHY simulates carbon assimilation and plant and soil respiration within a full energy and water balance. Although we prescribe leaf area index (LAI) to the model, this version of BETHY has an optional leaf area dynamics module for prognostic LAI as described in Knorr et al. (2010). BETHY represents variability in the physiology and ecology of plant classes by 13 PFTs (see Table 1) originally based on classifications by Wilson and Henderson-Sellers (1985). Each model grid cell may consist of up to three PFTs as defined by their grid cell fractional coverage.

Table 1PFTs defined in BETHY and their abbreviations.

Download Print Version | Download XLSX

SCOPE (version 1.53) is a vertically integrated (one-dimensional, 1-D) radiative transfer and energy balance model with modules for photosynthesis and chlorophyll fluorescence (van der Tol et al., 2009). It utilizes a canopy radiative transfer scheme based on the Scattering by Arbitrarily Inclined Leaves (SAIL) model (Verhoef, 1984) and the leaf radiative transfer model of Fluspect (Miller et al., 2005), which is based upon the optical properties of leaves (Jacquemoud and Baret, 1990). A limitation of this SCOPE version is that it lacks a water balance and only accounts for vertical variation in canopy properties, not horizontal variation. We note that a recent update has included a water balance and water stress in SCOPE, although this was only tested at a semiarid grassland site (Bayat et al., 2019).

While van der Tol et al. (2009, 2014) provide a more comprehensive description of the SCOPE model, we provide a brief description of the link between SIF and GPP. During the iterative calculation of the thermal radiative transfer and energy balance modules, the photosynthesis and chlorophyll fluorescence quantum efficiency (ϕ_F) of each canopy element are calculated. This includes the leaf biochemistry module, which simulates the photosynthetic rate as the minimum of two potentially limiting reaction rates (see Collatz et al., 1991 for C₃ plants and Collatz et al., 1992 for C₄ plants). Inputs to the leaf biochemistry module include APAR, relative humidity, temperature, CO₂ concentration, O₂ concentration and leaf physiological parameters (e.g., carboxylation capacity). The ϕ_F for each canopy element is also calculated within this module. To determine the ϕ_F, the fate of absorbed quanta via PQ and NPQ must be determined. First, the photochemical yield (ϕ_P) is determined from the quantum requirement of the photosynthetic dark reactions, i.e., the electron transport rate (J_e), where J_e is calculated as

where A_g is the gross photosynthetic rate (i.e., excluding dark respiration), C_i is the intercellular CO₂ partial pressure and Γ^* is the CO₂ compensation point, and effcon is a variable based on the electron requirements and assumptions on the processes limiting electron transport (in SCOPE, C₃ plant $\mathrm{effcon}=\mathrm{1}/\mathrm{5}$; C₄ plant $\mathrm{effcon}=\mathrm{1}/\mathrm{6}$). Photochemical yield is determined by the ratio of J_e to the total absorbed flux of electrons (J_APAR). The ϕ_P is calculated, based on the Genty et al. (1989) relationship (see van der Tol et al., 2014), as follows:

where J_APAR is the flux of electrons absorbed by photosystem II, assumed to equal to half the APAR. In this study, the APAR driving biochemistry is only that absorbed by chlorophyll, which differs from the original SCOPE model but is consistent with understanding of light-harvesting (Fleming et al., 2012). Equation (2) imposes the condition that the flux of electrons produced from photochemistry must equal those consumed by the photosynthetic dark reactions (van der Tol et al., 2014). The remaining quanta are distributed between chlorophyll fluorescence and NPQ, where NPQ is further split into constitutive thermal dissipation (constitutive NPQ) and energy-dependent, regulated thermal dissipation (regulated NPQ).

From the Genty et al. (1989) relationship, the steady-state ϕ_P equals the ratio of variable fluorescence to total fluorescence: ${\mathit{\varphi }}_{\mathrm{P}}=(F_{\mathrm{m}}^{\prime }-F_{\mathrm{t}})/F_{\mathrm{m}}^{\prime }$, where F_t is the steady-state fluorescence and $F_{\mathrm{m}}^{\prime }$ is the maximal fluorescence under a saturating pulse, indicating regulated NPQ. This evolved from decades of research using pulse amplitude modulation fluorescence measurements and theory (Baker, 2008). This relationship can be rearranged to the following:

Therefore, to obtain ϕ_F, a formulation for $F_{\mathrm{m}}^{\prime }$ is required. Understanding of the mechanisms driving $F_{\mathrm{m}}^{\prime }$ are not yet sufficient for a process-based model applicable in a canopy-scale steady-state photosynthesis model; however, van der Tol et al. (2014) showed that its variability could be captured using an empirical formulation. Here, we use the empirical fit to the drought data (Flexas et al., 2002; van der Tol et al., 2014). It is known that regulated NPQ is controlled by biochemical feedbacks (Zaks et al., 2013). It is therefore calculated using an empirically derived equation and a variable that describes the strength of the feedback termed the relative light saturation of photosynthesis, defined as $\mathrm{1}-{\mathit{\varphi }}_{\mathrm{P}}/{\mathit{\varphi }}_{\mathrm{P}}^{\mathrm{0}}$ (see van der Tol et al., 2014), where ${\mathit{\varphi }}_{\mathrm{P}}^{\mathrm{0}}$ is the maximum potential photochemical yield with typical values of 0.83 (Björkman and Demmig, 1987). Constitutive NPQ is also calculated, but it is known to be low and relatively constant, although the model does include a high temperature correction (van der Tol et al., 2014). Chlorophyll fluorescence quantum yield can thus be calculated by Eq. (3).

The photosystem I (PSI) and photosystem II (PSII) fluorescence spectra are calculated by the Fluspect module based upon the canopy structure, irradiance, biophysical properties (i.e., leaf composition; pigment concentrations, mesophyll structure and senescent and dry matter contents) and fluorescence quantum efficiency values for low-light unstressed conditions. Only the PSII fluorescence spectra is adjusted for regulatory feedbacks, as PSI fluorescence is considered to be relatively low and constant (Porcar-Castell et al., 2014). This is modeled by scaling the PSII spectra with the ratio ϕ_F from Eq. (3) (sometimes denoted by η_II) to the low-light unstressed quantum yield (sometimes denoted by η_II(0)). Re-absorption and scattering of fluorescence within the canopy is calculated by a separate routine (van der Tol et al., 2009). This is wavelength-dependent and occurs based on leaf composition and canopy structure. The fluorescence of canopy elements are then numerically integrated over canopy depth and orientation to determine the top-of-canopy SIF, similarly performed for leaf photosynthetic rates to determine GPP.

Overall, the modeled link between SIF and GPP occurs via the above equations. Therefore, variables (e.g., input parameters, environmental variables) that affect the photosynthetic rate will also affect SIF via ϕ_F. This includes variables that affect APAR, as APAR is an input to the leaf biochemistry module. However, APAR not only modulates ϕ_F, but has the additional, perhaps more significant effect of scaling the fluorescence spectra. Furthermore, variables such as leaf composition or canopy structure can influence the escape probability of fluorescence emission by re-absorption and scattering.

Between BETHY-SCOPE v1.0 used in Norton et al. (2018) and v1.1 used here, the key changes are (i) the correction of an error in the Fluspect module where the fluorescence quantum efficiency (f_qe) for PSI and PSII were set to be equal, while SCOPE v1.53 sets f_qe for PSI to be one-fifth that of PSII, and (ii) the leaf biochemistry module is now driven by green APAR (as mentioned above), rather than total APAR that is used in SCOPE v.153. Overall, in BETHY-SCOPE the canopy radiative transfer, energy balance and leaf biochemistry schemes of BETHY have been replaced by the corresponding schemes in SCOPE. The spatial distribution, vegetation (PFT) characteristics and carbon balance are handled by BETHY. SCOPE therefore takes climate forcing (meteorological and radiation data) and spatial information from BETHY and returns GPP, enabling process-based global simulations of GPP and SIF.

2.2 BETHY-SCOPE parameters

In this data assimilation system, the quantities to be optimized are the biophysical parameters that relate to SIF and GPP (see Table A1 in Appendix A). Parameters can be either global or spatially differentiated by PFT. PFT-dependent parameters enable differentiation between biophysical traits. Two key parameters for this study, the maximum carboxylation rate at 25 ^∘C (V_cmax) and chlorophyll a∕b content (C_ab), are considered PFT-dependent. The V_cmax parameter is used in most process-based terrestrial biosphere models as it is a parameter of the photosynthesis model of Farquhar et al. (1980). The C_ab parameter is a parameter specific to the SCOPE model and an important component of the canopy radiative transfer scheme as it strongly influences both SIF and APAR.

In total there are 41 parameters that are optimized by the data assimilation system. The uncertainty associated with each of these parameters is represented by a Gaussian probability density function (PDF). The mean and standard deviation for the prior parameters are shown in Table A1. Choice of the prior mean and uncertainty follow those used in previous studies (Kaminski et al., 2012; Knorr et al., 2010; Koffi et al., 2015). For new parameters that are not well characterized (e.g., SCOPE parameters) we assign relatively large prior uncertainties and mean values in line with the default SCOPE parameters and with Koffi et al. (2015) and Norton et al. (2018). An exception is the C_ab parameters, which are assigned higher prior values than Norton et al. (2018), more in line with physiological understanding.

Parameters exposed to the data assimilation system are chosen based on previous sensitivity tests such as those performed by Verrelst et al. (2015) and Norton et al. (2018). This includes leaf composition parameters such as C_ab, leaf dry matter content (C_dm) and leaf senescent material fraction (C_s). Also included are structural parameters such as leaf distribution function parameters (LIDFa, LIDFb), vegetation height (hc) and leaf mesophyll structure, the prior values for these were obtained from literature values and are assigned to groups of PFTs that we assume have a generally similar structural form (see Table A1). Physiological parameters are also incorporated, including V_cmax and Michaelis-–Menten constants of Rubisco for CO₂ (K_C) and O₂ (K_O). Additionally, the photosynthetic kinetic parameter for the maximum oxygenation rate (V_omax) is included. Given the uncertainty of V_omax and its importance for modeling GPP (von Caemmerer, 2000), this may be an important parameter to consider and is given by its ratio with V_cmax, $a_{V_{\mathrm{o}},V_{\mathrm{c}}}$. Given this parameter also affects the relative specificity of Rubisco (S_c∕o) we calculate S_c∕o explicitly following von Caemmerer (2000) which differs from the original SCOPE model.

2.3 Satellite SIF observations

We use satellite SIF observations from the NASA Orbiting Carbon Observatory-2 (OCO-2) (Sun et al., 2018). Launched in July 2014, OCO-2 operates in a sun-synchronous orbit with an overpass at approximately 13:30 local time and a repeat cycle of 16 d. Collecting approximately 24 spectra per second, it has relatively high data density within the field of view. OCO-2 has a ground-pixel size of 1.3×2.25 km² and a total swath width of 10.6 km. Full spatial mapping of SIF is therefore not possible with OCO-2. However, the high spectral resolution of OCO-2 allows for robust and accurate SIF retrievals (Frankenberg et al., 2014; Sun et al., 2018).

Alternative satellite SIF datasets are also available, including from the GOME-2 and GOSAT instruments (Frankenberg et al., 2011a; Guanter et al., 2012; Joiner et al., 2011). There are benefits and pitfalls in using these alternative data. For example, GOME-2 and GOSAT provide longer time series going back to 2007 and 2009, respectively. GOME-2 also provides better spatial mapping compared with OCO-2. However, there are known issues of sensor degradation with GOME-2 (Zhang et al., 2018). The advantage of the OCO-2 satellite is that it collects 8 times more spectra and has a higher spectral resolution providing more robust and data dense observations (Frankenberg et al., 2014; Sun et al., 2018). We note that a formal comparison of these other datasets is outside the scope of this study, but a recent comparison of TROPOMI and OCO-2 showed strong agreement (Köhler et al., 2018).

We use the OCO-2 processed SIF-lite data files. For details on the retrieval algorithm for the SIF data see Frankenberg et al. (2014) and Sun et al. (2018). These data are gridded at 2^∘ × 2^∘ spatial resolution, equivalent to the model grid resolution. We exclude soundings collected over water as determined by the corresponding International Geosphere-Biosphere Programme (IGBP) land classification index (Friedl et al., 2010). We use instantaneous SIF at 757 nm and only soundings taken in nadir mode. Data are also available at 771 nm; however, the signal at 757 nm is stronger (Sun et al., 2018), and thus we only consider that signal. The annual mean OCO-2 SIF for 2015 is shown in Fig. 1.

There are potential limitations in using OCO-2 for global mapping due to spatial coverage of the observations and sampling bias of biomes, particularly if it differs from the assumed biome types used in the model. We assessed the similarity between the sampled IGBP land classification index of the OCO-2 soundings and the BETHY-SCOPE PFTs to evaluate this limitation. Considering the differences in vegetation classifications this is a qualitative test. Qualitatively, the occurrence of IGBP biome sampling appears to be similar to the BETHY-SCOPE PFTs. Moreover, Frankenberg et al. (2014) showed that despite the limited spatial coverage of OCO-2 it provides a representative sampling of 1^∘ × 1^∘ grid cell averages. We therefore do not perform any further filtering of the data. Future studies may benefit from evaluating this in more detail and performing the assimilation using observations at a PFT-specific level rather than the grid cell level as is applied here.

2.4 Data assimilation system

To assimilate SIF into BETHY-SCOPE we require a minimization algorithm, cost function and error propagation method. A variety of techniques are available for the optimization of terrestrial biosphere models and reviews are available (Fox et al., 2009; Kaminski et al., 2013; MacBean et al., 2016; Trudinger et al., 2007).

We utilize a probabilistic framework whereby quantities (e.g., observations, model state variables, model process parameters) are represented by their PDF. These quantities are treated as Gaussian, and thus can be described by their mean and standard deviation. For the model parameters the mean is denoted by x and error covariance matrix by C_x. We denote the prior parameter vector and covariance matrix by x₀ and ${\mathbf{C}}_{x_{\mathrm{0}}}$, respectively, and the posterior parameter vector and covariance matrix by x_post and ${\mathbf{C}}_{x_{\mathrm{post}}}$, respectively. For the observations the mean is denoted by d. The error covariance matrix in observation space, denoted by C_d, combines errors in the observations and in their simulated counterpart, i.e., model (Kuppel et al., 2013). Quantification of model error can be performed through an assessment of model–observation residuals following optimization (e.g., Kuppel et al., 2013). We assess potential model errors in this study; however, we do not explicitly account for this error in the propagation of errors onto GPP; hence C_d accounts only for errors in the observations. We point out that the uncertainty is embodied in the error covariance matrices and that diagonal elements represent the variance of the quantities while off-diagonal elements represent error covariances between quantities.

2.5 Experimental setup

In this study BETHY-SCOPE is run for the year 2015. This constitutes the optimization (or calibration) period. We then assess the optimized model performance against independent OCO-2 observations outside of the optimization period from September to December 2014.

The model is run on a 2^∘ × 2^∘ grid resolution. Climate forcing data are provided in the form of monthly meteorology (precipitation, minimum and maximum temperatures, and incoming solar radiation) obtained from the WATCH/ERA Interim dataset (WFDEI; Weedon et al., 2014). These are used to derive average diurnal cycles of meteorological forcing. Thus, a single average diurnal cycle of meteorological forcing for each month is used to simulate photosynthesis and fluorescence. This allows the computation of SIF at the equivalent overpass time as the satellite data (13:00–14:00 local time). The model SIF is also calculated at the equivalent wavelength to the OCO-2 SIF data (757 nm). Atmospheric CO₂ concentration is set to the 2015 annual average. LAI is prescribed to the model using the MODIS improved LAI dataset (Yuan et al., 2011). The LAI is averaged at the model 2^∘ × 2^∘ grid resolution and for each grid cell it is split between PFTs using the model PFT grid cell fractional coverage.

2.6 Global GPP products for comparison

To assess the SIF-optimized global GPP we compare the model prior and posterior GPP to other global GPP products. The first dataset for comparison is an upscaled product based on site-level measurements termed FLUXCOM GPP (Tramontana et al., 2016). The FLUXCOM GPP product uses various machine-learning techniques to empirically upscale flux tower data using remote sensing and meteorological data as the predictor variables. Here, we use the ensemble average of the FLUXCOM GPP product that uses remote sensing data exclusively (Tramontana et al., 2016). The second dataset for comparison is an ensemble of 11 global dynamic vegetation models forced with equivalent climate fields and atmospheric CO₂ concentration that were used to investigate trends in sources and sinks of CO₂ (TRENDY; Sitch et al., 2015). It is important to note that global GPP is highly uncertain and that both the FLUXCOM GPP and TRENDY GPP estimates are based on their own model assumptions and/or sparse measured data (Anav et al., 2015). Therefore, these data are used to evaluate whether the SIF assimilation results in global patterns of GPP that align with the current understanding and not for extensive validation purposes.

3.1 Assimilation with SIF

Here we show how the model compares with the observed SIF (shown in Fig. 1) for the prior and posterior cases and for the calibration and validation periods. The goodness of fit between modeled and observed SIF is assessed using multiple metrics. The ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ fit is a key metric (see Eq. 7). Differences between the model and observations (“residual”) and the squared residual normalized by the observational variance (“mismatch”) are also used. The mismatch is a measure of the difference between the model and observations accounting for observational uncertainties, indicating the contribution grid cells make toward the cost function. The model fit during the calibration period is presented in more detail as there is more data. The model fit during the validation period provides a more stringent test of the assimilation performance. We then show the performance of an additional model simulation testing seasonal variation in parameters.

Figure 1Annual mean observed SIF from the OCO-2 satellite for 2015.

3.1.1 Calibration

The model prior SIF (SIF_prior) over the calibration period yields a global ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ fit of 2.45. Large residuals in the annual mean, shown in Fig. 2, are evident across the globe, ranging from −0.87 to +1.20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$. Generally, SIF_prior overestimates observed SIF across regions dominated by tropical forest (e.g., the Amazon, western equatorial Africa and Maritime Continent), boreal forest (parts of North America and Eurasia) and semiarid regions (e.g., central Australia, central Asia and southern Africa). SIF_prior tends to underestimate observed SIF across the rest of the land, in particular for regions dominated by croplands (e.g., American Midwest, parts of Europe, India, eastern Asia), mixed forests (across Europe and Asia), and grassland and savanna regions (e.g., the African savanna and Brazilian Highlands of South America). Latitudinal averages, shown in Fig. 4, also show these spatial patterns for SIF_prior. Overestimation of observed SIF is seen over the central tropics between 15^∘ S and 5^∘ N, a region dominated by tropical evergreen forest (TrEv), whereas there is significant underestimation of observed SIF over the Northern Hemisphere, particularly during northern summer (see Fig. 5).

Figure 2Annual mean residual between model SIF_prior and observed SIF for 2015.

Figure 3Annual mean residual between model SIF_post and observed SIF for 2015.

Following the assimilation, the model shows a considerably better fit to the calibration data. The global ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ fit is strongly reduced from 2.45 to 1.01. This is close to the optimal value of 1, demonstrating the ability of the optimized model to fit the observed patterns of SIF and validating our chosen uncertainties as far as is practicable, including the choice of the scaling used to calculate observational uncertainties in Eq. (4). Annual mean residuals between model posterior SIF (SIF_post) and the observations, shown in Fig. 3, range between −0.58 and +0.45 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, considerably smaller than SIF_prior. The spatial patterns of posterior residuals (Fig. 3) show that the regional patterns of residuals broadly exhibit the same sign as the prior case albeit with a significantly reduced magnitude. For the latitudinal averages, SIF_post is remarkably close to the observed SIF for the annual average (Fig. 4). However, discrepancies are evident for the northern summer averages for both SIF_prior and SIF_post (Fig. 5), where observed SIF is underestimated in the Northern Hemisphere and overestimated in the Southern Hemisphere.

Latitudinal sums of the mismatch between the model and the observations (bar charts in Figs. 4 and 5) also show a significant reduction (i.e., improvement in fit) following the assimilation. Comparison of the gray and green bars, indicating the respective prior and posterior mismatch, show that this improvement in fit occurs over all latitudes at the annual and northern summer timescales. The total annual mismatch between SIF_post and the observations is about 40 %–80 % smaller across the latitudes between 40^∘ S and 60^∘ N relative to SIF_prior.

Despite the strong improvement in fit, SIF_post tends to underestimate large observed SIF values >1.0 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, shown in the Fig. S4. A linear regression line between the observed SIF and SIF_post has a slope of 0.67. Furthermore, while the observed SIF for any given month and grid cell can reach up to 2.29 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, SIF_prior does not exceed 1.98 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ and SIF_post does not exceed 1.43 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$. The SIF_prior does not show this systematic underestimation, however, but it has a poorer global fit (Fig. S3). We note that these large observed SIF values occur mostly over the tropics and the northern midlatitudes during the peak growing season.

From Fig. 3 it appears that SIF_post overestimates observed SIF over arid regions (e.g., the Sahara, Atacama and Namib deserts, central Australia and central Asia). This occurs primarily because the observed SIF is slightly negative (see Fig. 1), potentially due to measurement noise or issues from the correction of constant error artifacts in the SIF retrieval (Sun et al., 2018). Negative SIF values are still considered in the assimilation system. However, they contribute little to the overall mismatch given the uncertainty in the SIF observations (see Figs. S8 and S9).

Figure 4Latitudinal averaged SIF and mismatch with observations. OCO-2 observations (black line), SIF_prior (gray dashed line, gray bars) and SIF_post (green line, green bars) for the annual average. Data are only shown for spatiotemporal points where OCO-2 observations are present.

Download

Figure 5Latitudinal averaged SIF and mismatch with observations. OCO-2 observations (black line), SIF_prior (gray dashed line, gray bars) and SIF_post (green line, green bars) for northern summer (June–August). Data are only shown for spatiotemporal points where OCO-2 observations are present.

Download

3.1.4 Fit to the seasonal cycle

We can also assess the seasonal cycle of SIF to determine how well the model simulates the amplitude of observed SIF. First, we assess how well the model replicates the seasonal amplitude of observed SIF across all spatial points. Second, we assess the seasonal patterns of SIF for a selection of case study regions in more detail. We avoid assessing the seasonal cycle of SIF aggregated at global or hemispheric scales as regional patterns of residuals can differ in sign and magnitude (e.g., see Fig. 3). The seasonal amplitude of observed SIF is calculated as the difference between the maximum and minimum SIF across the year for each grid point. To increase confidence that the observations really capture the seasonal cycle, we only assess grid points with at least 8 months of observed SIF data. In doing so, most regions north of 60^∘ N are excluded due to limited SIF observations. We do not assess the timing of the seasonal cycle (e.g., start and end of the growing season) as this is largely driven by LAI which is prescribed and therefore fixed in this study.

The comparison of the seasonal amplitude of observed SIF against SIF_prior, SIF_post and SIF_post,seas is shown in Fig. S10. The model underestimates the observed seasonal amplitude in all cases. With a perfect match to the observed seasonal amplitude the model would follow the 1:1 line and the slope of a linear regression line would equal 1. However, we find that the slope is 0.21 for SIF_prior, 0.35 for SIF_post and 0.43 for SIF_post,seas. Spatial points with the largest seasonal variations in observed SIF also exhibit the largest model-observed mismatch (Fig. 6).

Figure 6Seasonal amplitude of observed SIF versus the model-observed mismatch for the SIF optimized (SIF_post) model. Mismatch is defined as the squared residual normalized by the observational uncertainty (expressed as variances). A positive curvi-linear relationship exists such that spatiotemporal points with a large model-observed mismatch generally also exhibit a large observed seasonal amplitude.

Download

For a more detailed assessment of seasonal patterns we investigate three case study regions: (i) the tropical forest of mainland southeast Asia, (ii) croplands in North America and (iii) the north African savanna (see Figs. S11–S15 for details). These regions are selected as they represent quite different biome types, exhibit varied SIF patterns and have relatively large posterior model-observed mismatch.

The tropical evergreen forest of mainland southeast Asia exhibits a clear seasonal cycle in observed SIF. The monthly mean observed SIF averaged over this region varies from a minimum of ∼0.6 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ in March to a maximum of ∼1.4 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ in August (see Fig. S11). Both SIF_prior and SIF_post exhibit seasonal cycles that differ strongly from the observations, with little change in the shape of the seasonal cycle from the prior to posterior simulations. Model SIF shows a minimum in July and two maximums in February and November, following the seasonal evolution of LAI (see Fig. S11), which we reiterate is prescribed. This results in strong negative temporal correlations between observed SIF and model SIF as over this region.

Croplands in North America are heavily managed landscapes with highly productive vegetation as indicated by the large observed SIF values during the growing season. Even with the monthly averages used here, observed SIF can exceed 2.0 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$. As presented earlier, the model SIF cannot match the seasonal amplitude of observed SIF and subsequently underestimates the maximum monthly SIF averaged over this region by almost 40 % (0.45 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$) (see Fig. S13). The fit is improved in SIF_post,seas (${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\left({\text{SIF}}_{\text{post}}\right)=\mathrm{4.62}$; ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\left({\text{SIF}}_{\text{post,seas}}\right)=\mathrm{3.10}$). In both cases the timing of seasonal maximum and senescence is simulated quite well, while the onset of the growing season is predicted to be too early.

The north African savanna exhibits a strong seasonal cycle in observed SIF. This region is dominated by grasslands and open forest, with a seasonality closely following the seasonal variation in precipitation. Averaged over the region, observed SIF varies from 0.07 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ in January to 0.77 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{\mu }{\mathrm{m}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ in September (see Fig. S15). However, SIF_post exhibits a much smaller seasonality with variation from 0.24 to 0.48 across the year, only 34 % of the observed seasonal amplitude. Temporal correlations are quite strong, however, as model SIF also reaches its peak in September.

3.2 Estimated parameters

The prior and posterior parameter mean values and associated uncertainties are shown in Table A1. In this data assimilation system the number of observations far outweighs the number of unknowns. This means that there is a substantial amount of observational information available to constrain parameter values, and thus they can shift from their prior values considerably even if given a relatively tight prior uncertainty. We can be more confident in parameters that see large reductions in uncertainty. Conversely, parameters with little reduction in uncertainty following optimization should be accepted cautiously.

Posterior V_cmax estimates range from 16 to 130 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. The lowest posterior rates occur for the TmpEv, C4Gr, Tund and Wetl PFTs, all equal or below 30 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. The highest posterior rates occur for the Crop and C3Gr PFTs, both exceeding 100 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. Out of 13 PFTs, 9 see an increase in V_cmax. Increases greater than 2 standard deviations occur for the PFTs TmpDec, EvCn, C3Gr, C4Gr and Tund, many of which dominate temperate regions. Latitudinal averages of V_cmax (Fig. S19) show a strong increase in temperate-zone V_cmax following the SIF assimilation, shifting it higher than V_cmax in the tropics. Zonally, the lowest V_cmax occurs in the cold-climate high latitudes. Most of these parameters see moderately strong uncertainty reductions (>30 %), indicating strong constraints from SIF.

Posterior C_ab estimates range from 8 to 38 µg cm⁻². Out of 13 PFTs, 9 see a decrease in C_ab. The highest posterior C_ab values occur for the TrEv, TmpEv, TmpDec, DecCn and Crop PFTs, all >25 µg cm⁻², while the lowest values occur for the EvShr, Tund and Wetl PFTs, all <15 µg cm⁻². Uncertainty reduction ranges from weak to moderately strong for the C_ab parameters, with a maximum constraint of 34 %. Latitudinal averages of posterior C_ab show a relatively low variance across different zones (Fig. S20), with posterior values being similar between the tropics and temperate midlatitudes, albeit with a distinct dip in drier subtropics and high latitudes. The leaf composition parameter for dry matter content, C_dm, increases by about 70 % and an uncertainty reduction of about 20 %. The leaf composition parameter for senescent material remains almost unchanged and shows only a very weak uncertainty reduction of about 1 %.

Parameters that control canopy structure and the leaf angle distribution see large deviations from their prior values. Some leaf angle distribution parameters, LIDFa and LIDFb, shift considerably. SIF is particularly sensitive to the LIDFa and LIDFb parameters, although this depends on which group of PFTs they pertain to, and so these parameters have uncertainty reductions of up to 90 %. Vegetation height for both trees and shrubs remain relatively unchanged, while vegetation height for grasses and crops see a decrease of over 4 standard deviations. However, the uncertainty reduction of the vegetation height parameters is very weak (<1 %). Despite these changes GPP is relatively insensitive to the canopy structure parameters.

3.3 Estimated GPP

The spatial patterns of posterior GPP and the changes following the SIF assimilation are shown in Figs. 7–11. Following the assimilation of SIF global GPP increases by about 39 Pg C yr⁻¹, from 127.6 to 166.7 Pg C yr⁻¹. This change shifts BETHY-SCOPE further from the TRENDY mean (142.4 Pg C yr⁻¹) and FLUXCOM GPP (103.3 Pg C yr⁻¹) estimates. The parametric uncertainty in global GPP is reduced by 38 % from ±7.4 to ±4.6 Pg C yr⁻¹ by the SIF assimilation.

Figure 7Spatial patterns of BETHY-SCOPE posterior annual mean GPP (GPP_post) for 2015.

Figure 8Change in annual mean GPP rate for 2015 following optimization with SIF relative to GPP_prior.

Spatially, increases in annual GPP are seen across much of the land surface as shown in Fig. 8 (see Fig. S17 for the percentage change in GPP). Declines in GPP are seen in dry tropical forests in parts of South America, Africa and mainland Asia (TrDec PFT; also shown in Fig. 9). Globally, TrDec-dominated forests see a decline in total GPP of 23 %. Two other biomes, TmpEv and DecShr, show a decline in global GPP; however, their contribution to global GPP is relatively small (see Fig. 9). The global GPP of TrEv forest biomes, including the Amazon, equatorial Africa and Maritime Continent, increase by about 20 %. In relative terms, the largest increases in global GPP occur for TmpDec, EvCn, C3Gr and C4Gr, all >60 %. We note that these changes in model GPP can differ in sign and magnitude from the changes in model SIF (see Figs. 8 and S18). This can occur as SIF and GPP have differing sensitivities to the underlying parameters, a result of the process-based approach. For example, the wet tropical forests (TrEv) and cold-climate conifer forests (EvCn) see an increase in GPP but decline in SIF.

Figure 9Annual GPP for 2015 per biome. Biomes are defined by aggregating model grid cells that have the same spatially dominant PFT as shown in the Fig. A1 in Appendix A.

Download

Figure 10Annual latitudinal averages of GPP_prior (gray line), GPP_post (green line), FLUXCOM GPP (orange line), TRENDY model average (light blue line) and TRENDY model spread given by the 10th and 90th percentiles (light blue shading).

Download

Averaged over latitudinal bands (Figs. 10 and 11) we can see the distinctive peak in GPP across the tropics at the annual timescale and secondary peak in GPP across the northern midlatitudes during the northern summer. With the assimilation of SIF the GPP across all latitudes increases. Some regions and seasons see larger changes. Across the central tropics (15^∘ S–5^∘ N; dominated by the PFT TrEv), GPP increases substantially. For this region, the prior and posterior estimates are near the high end of other estimates, with the posterior GPP exceeding the 90th percentile of TRENDY models (blue shading). FLUXCOM GPP is very low in the central tropics, but we note that this product is not expected to be representative of the tropics given the sparsity of the flux tower network there (Tramontana et al., 2016). The northern extratropics (30–60^∘ N) show a general increase in BETHY-SCOPE GPP, with the SIF assimilation shifting it to the higher end of other estimates. This is particularly strong during northern summer (Fig. 11). While the prior GPP in this region is within the TRENDY model range and close to the FLUXCOM GPP, the SIF assimilation results in a posterior that exceeds the 90th percentile of the TRENDY model range. North of 65^∘ N the prior closely matches FLUXCOM GPP, but the posterior sees an increase which brings it more in line with the TRENDY model average. There is also a distinct difference between the FLUXCOM GPP product and all models north of 75^∘ N, with FLUXCOM GPP being higher. BETHY-SCOPE posterior GPP over the southern latitudes south of 15^∘ S is generally within the TRENDY model range. In this region, the prior GPP is near the bottom of the TRENDY model range, although this shows closer similarity to the FLUXCOM GPP.

Figure 11Northern summer (June–August) latitudinal averages of GPP_prior (gray line), GPP_post (green line), FLUXCOM GPP (orange line), TRENDY model average (light blue line) and TRENDY model spread given by the 10th and 90th percentiles (light blue shading).

Download

Figure 12Per PFT LUE_GPP for the prior and posterior simulations. Values were determined by monthly average GPP divided by monthly averaged APAR, then averaged to annual scales thus averaging over seasonal variations. Note that the theoretical maximum is 0.08 mol C mol photons⁻¹ (Waring et al., 2016).

Download

A useful metric for patterns of global productivity is the ratio of GPP between different regions. These ratios are summarized in Table B1 in Appendix B. The ratio of the tropics (30^∘ S–30^∘ N) to the extratropics (south of 30^∘ S and north of 30^∘ N) declines following the SIF assimilation, due to a relatively larger increase in extratropical GPP compared to tropical GPP. This shifts the ratio of tropical to extratropical GPP from a prior of 2.47 to a posterior of 1.94, which is substantially closer to patterns of the FLUXCOM (1.90) and TRENDY (1.93) mean. Similarly, the ratios of the tropics to the boreal region (north of 55^∘ N), tropics to the temperate region (south of 30^∘ S and north of 30–55^∘ N) and the temperate to boreal region converge toward the FLUXCOM values (see Table B1).

We also note an improvement in the correlation between the BETHY-SCOPE estimate and the FLUXCOM GPP over North America (see Figs. B1 and B2 in Appendix B) and Europe (data not shown), two regions where FLUXCOM GPP has considerably more training data. Over North America the correlation improves from a prior R²=0.80 to a posterior of R²=0.89, while over Europe the correlation improves from a prior R²=0.76 to a posterior of R²=0.86. Despite this improvement in match between the patterns, the posterior slope is 1.6 for North America and 1.8 for Europe, indicating that the magnitude of posterior monthly GPP is larger than that of FLUXCOM GPP.

Changes in GPP due to changes in parameter values can be broken down into changes in intercepted radiation (APAR) and canopy photosynthetic light-use efficiency (LUE_GPP). The LUE_GPP is calculated as the annual average of the ratio between monthly GPP to monthly APAR. Overall, the majority of biomes see an increase in LUE_GPP following the SIF assimilation, shown in Fig. 12. Just three biomes, TrDec, TmpEv and DecShr, see a decline in LUE_GPP. These biomes are the same that see a decline in annual GPP (comparing Figs. 8 and B4). This can be related back to the general increase in V_cmax for most PFTs (Table A1) and latitudes (Fig. S19). Changes in APAR are smaller in relative terms and show distinct regional differences (Fig. B6). With the exception of low-productivity arid regions, the largest percentage change in APAR occurs for the high-latitude tundra biome with approximately a 20 % increase in APAR, due primarily to an increase in C_ab. Wet tropical forests see a decline in APAR of about 5 %, while drier tropical biomes (e.g., Brazilian Highlands, north African savanna) see an increase of <5 %. Other regions show only minor shifts in APAR.