Inferring CO₂ fertilization effect based on global monitoring land-atmosphere exchange with a theoretical model

We used a 104 flux tower sites (770 site years) (Ichii et al 2017, Pastorello et al 2017) covering the period of 1990–2014 across a number of biome types (forests, grasslands, savanna, shrub, wetlands, tundra, and cropland) (figure S1, table S1); 66 sites were forests, whereas 38 sites were non-forested. The data consisted of two data sets, one for Asia, the other covering the rest of the globe. The data set for Asia was an extended version of our previous data set (Ichii et al 2017), which consisted of 60 EC observation sites and these were derived from five databases, namely AsiaFlux (http://asiaflux.net), European Flux Database (http://europe-fluxdata.eu), Forestry and Forest Products Research Institute (FFPRI; http://2.ffpri.affrc.go.jp/labs/flux/), National Institute Agro-Environmental Studies (NIAES), and Arctic Observatory Network (AON; http://aon.iab.uaf.edu/). The data set for the rest of the globe was 44 sites from the FLUXNET2015 (http://fluxnet.fluxdata.org/data/fluxnet2015-dataset/; Pastorello et al 2017), where we selected 39 sites that had more than ten years of data and an additional five sites with less than ten years of data to capture South American and African regions.

We used the original EC measurements of NEE that were not post-processed by principal investigators or their respective Asian databases, and partitioned GPP from FLUXNET2015. We post-processed the Asian data for all sites using a standard protocol for data quality control, and flux partitioning (Ichii et al 2017). Quality control was conducted using a spike detection method (Papale et al 2006) and a filtering with the friction velocity threshold (Reichstein et al 2005). GPP was estimated as the difference between net ecosystem exchange (NEE) and ecosystem respiration (RE). Daytime RE was based on exponential relationships (Lloyd and Taylor 1994), which were determined for each day using nighttime data with a 29-day moving window. GPP and RE were determined as the mean of 100 values obtained from an ordinary bootstrapping procedure. Further details are available in A3 in the supplemental material. For FLUXNET2015 data, we used partitioned GPP based on nighttime data (GPP_NT_VUT_REF) as this approach was similar to the processing for the Asia data set. No gap-filled flux data were used for model optimization in both data sets. The current study focused on C3 photosynthesis, because C4-dominated ecosystems were limited in our data set and a much smaller CO₂ fertilization effect is expected in C4 ecosystems than C3-dominated ecosystems (Collatz et al 1992).

We used the canopy photosynthesis model (iBLM-EC version 2; further for details can be found in the Supplemental Material) to quantify the CO₂ fertilization effect (Ueyama et al 2016, 2018). The model contains coupled effects among the following sub-models: biochemical photosynthesis (Farquhar et al 1980), stomatal conductance (Ball et al 1987), aerodynamic conductance (Lhomme et al 2000), sun/shade radiative transfer including clumping effect (Ryu et al 2011), and vertical nitrogen transfer within the canopy (Lloyd et al 2010). In the biochemical photosynthesis model, photosynthesis is modeled as the minimum rate of Rubisco-limited photosynthesis and RuBP-limited photosynthesis. Rubisco-limited photosynthesis is sensitive to intercellular CO₂ concentration, whereas RuBP-limited photosynthesis is sensitive to light. Stomatal conductance is modeled based on a semi-empirical relationship to photosynthetic rate, CO₂ concentration at the leaf surface, and humidity (Ball et al 1987). Both photosynthesis and stomatal conductance are further calculated separately for sunlit and shaded leaves, where sunlit and shaded portions are calculated based on a radiative transfer model (Ryu et al 2011) that uses leaf area index (LAI) and clumping index as inputs (He et al 2012, Pisek et al 2015).

Input data for the model include half-hourly or hourly meteorology (wind speed, photosynthetically photon flux density, air temperature, relative humidity, rainfall, atmospheric pressure, net radiation, and ground heat flux), turbulent fluxes (friction velocity, sensible and latent heat fluxes, and GPP), [CO₂], LAI, and growth temperature (defined as a mean air temperature over the preceding thirty days). We did not apply the energy imbalance correction for inputs of sensible and latent heat fluxes except for a sensitivity analysis, because uncertainties in this assumption were negligible as shown in Supplemental Material. We rejected data due to wet conditions during rain and within one hour after rain, because uncertainties in this case were larger of the acceptance threshold 10% of estimated CO₂ fertilization effect, as shown in supplemental material. Transpiration was partitioned from measured latent heat flux by subtracting evaporation. The evaporation was estimated from potential evaporation at the soil surface for forests (Ryu et al 2011) solving net radiation at forest floor. The evaporation was estimated using LAI by known ratios between ET and transpiration for grassland, tundra, cropland (Wang et al 2014) and rice paddy (Sakuratani and Horie 1985). The direct and diffuse portions of radiation were partitioned based on the method (Weiss and Norman 1985) for solving the sun/shade radiative transfer model. To avoid the effects of inconsistent calibration, arising from inaccuracies in [CO₂] measurements at the observation sites, and for consistency among the optimization and CO₂ fertilization experiment, [CO₂] was derived from CarbonTracker 2015 (CT2015) (Peters et al 2007), with selected pixels centered on each site's geographic position. Since [CO₂] by CT2015 was available since 2000, [CO₂] prior to 2000 was estimated by subtracting the [CO₂] difference between values from 2000 and values from a target year at the Mauna Loa observatory from the CT2015 data for 2000. LAI in each site was from field observations if available (table S1). If time-variant LAI was not available, LAI was derived from MODIS MOD15A2H after smoothing the signal using TIMESAT (Eklundh and Jönsson 2015) and adjusting the mean annual maximum using site observations from literature (table S1).

We derived the parameters of the sun/shade canopy photosynthesis model from the constraint optimization (table S3) using the observed data to ensure that model accurately predicted the responses of EC-derived GPP, canopy-integrated gs, and thus transpiration to rising [CO₂]. Parameters that are sensitive to CO₂ fertilization, i.e. Vc_max25 and J_max25 in the biochemical photosynthesis model (Farquhar et al 1980) as well as m_bb and b_bb in the stomatal conductance model (Ball et al 1987), were determined using EC-based GPP and transpiration data. The parameters were determined based on a global optimization method (SCE-UA; Duan et al 1992) for each day and each site with an eight-day moving window (figure S3). We only used successfully converged parameters for estimating the CO₂ fertilization effect. We calculated mean and standard deviations of the parameters based on ten optimization runs from randomly generated initial values. We rejected the data when standard deviation of determined parameters was greater 10% of their absolute value for minimizing equifinality in parameterization (Medlyn et al 2005). Thus, uncertainties of the parameters associated with the optimization were less than 10% of the value. For quantifying the range of uncertainties from the biochemical photosynthesis model and parameterizations that were not directly constrained using observations, we quantified the CO₂ fertilization effect using six different submodels (Collatz et al 1991, de Pury and Farquhar 1997, Bernacchi et al 2001, 2003, Kosugi et al 2003, Kattge and Knorr 2007, von Caemmerer et al 2009). We used the ensemble mean and standard deviations of the estimated six CO₂ fertilization effects, which yielded a coefficient of variation of 10.3% across the six different parameterizations of the biochemical model (figure S4). Root mean square errors in half-hourly GPP and latent heat flux were 2.60 μmol m⁻² s⁻¹ and 34.6 W m⁻², respectively; slopes and R² between observations and the parameterized model also supported the validity of the optimization (figure S5).

Since the actual parameters for CO₂ fertilization processes (CO₂ demand and supply functions for photosynthesis; von Caemmerer et al 2009) were constrained each day through the optimization, the parameterized model can predict a theoretical sensitivity to rising atmospheric [CO₂] for a given day. The CO₂ fertilization effect was calculated at the half-hourly or hourly time step, and then averaged on an annual basis using daytime data when photosynthetically photon flux density (PPFD) was greater than 500 μmol m⁻² s⁻¹. The quantified CO₂ fertilization effect represents the changes in GPP and gs by changes in [CO₂] whilst allowing environmental conditions (e.g. climate and LAI) to change over time (Ueyama et al 2016). Here, we used [CO₂] in the year 2000 as a baseline (369.55 ppm of the annual CO₂ concentration at Mauna Loa) (Le Quéré et al 2018) and evaluate the effects of changing [CO₂] before and after 2000.

We separated two processes related to the changes in GPP and gs, namely the passive response by a hyperbolic response of photosynthesis to [CO₂] and the ecophysiological acclimation, known as down regulation. We defined the ecophysiological acclimation as the interannual/decadal variations in the ecophysiological parameters throughout observation periods. Based on this definition, we estimated the contributions of the ecophysiological acclimation to changes in GPP and gs by isolating the components associated with changes in interannual/decadal variations in the ecophysiological parameters (Vc_max25, J_max25, m_bb and b_bb). To examine the possible change associated with the ecophysiological acclimation, another experiment was conducted for sites having more than four years of data, where the CO₂ fertilization experiment was conducted using median seasonality in the parameters instead of daily parameters in each year. CO₂ fertilization effects between the two experiments would be different, if interannual variations or a directional shift in the parameters played an important role in determining the magnitude of the fertilization effect. The comparison was done for the years after 2000 when accurate ancillary inputs, such as LAI and [CO₂], were available.

We tested how the EC data constrained the CO₂ fertilization effect using data from 14 sites selected from the FLUXNET2015 data set (table S2). As a non-constrained experiment, we input biased parameters (Vc_max25, J_max25, and m_bb) into the model. We added or subtracted one standard deviation of the parameter distribution within a plant functional type to the determined parameters, and then estimated the range of uncertainties for the non-constrained simulation. In the non-constrained simulation, relative uncertainties of the CO₂ fertilization for GPP, gs, and iWUE were 23% ± 11%, 60% ± 58%, and 4% ± 5%, respectively. This indicates that the EC data effectively reduced the uncertainties in the CO₂ fertilization for GPP and gs.

For evaluating the CO₂ fertilization effect at the global scale, the effect examined at the site level was upscaled using a random forest regression, which is a machine-learning based regression. Changes in fluxes per change in [CO₂] (unit of % per ppm) were estimated for each site that contained more than four years of data (66 sites; table S1). The changes in fluxes were then modeled using a random forest regression by growing season and annual mean air temperatures, annual sum of precipitation, seasonal maximum of LAI, land cover (the alternative of forest or non-forest), and growing season length defined as number of days when mean air temperature was greater than 5 °C. We applied a 5 × 2 nested cross-validation with random search to tune generalized hyper parameters of the random forest regression; the nested cross validation effectively splits train, validation, and test sets for generalization. We measured the generalization error using the nested cross validation rather than using hold-out data due to the limited available data. For estimating possible uncertainties in the model tuning, we constructed 20 random forest regressions which had different hyper parameters using the 5 × 2 cross-validation scheme from 20 different initial parameters. The estimated R² score was 0.61 ± 0.05 for the effects for GPP and 0.78 ± 0.05 for gs.

The CO₂ fertilization effect was upscaled to the globe using the random forest regressions and gridded climate and satellite remote sensing data. For the global analyses, we calculated the CO₂ fertilization effect at an annual timescale for each grid that had a 0.25° × 0.25° spatial resolution. The change in GPP and ET associated with rising [CO₂] was estimated for the globe from 2000 to 2014 at the annual timescale. The mean climatology of annual mean air temperature and annual sum of precipitation from 2000 to 2014 was created based on Ichii et al (2013) by merging CRU-TS climate data and National Centers for Environmental Prediction (NCEP)/the National Center for Atmospheric Research (NCAR) reanalysis (NCEP/NCAR reanalysis) data with the quarter-degree spatial resolution (Ichii et al 2013). The maximum LAI was derived from MODIS standard product (MOD15A2H; collection 6) at the quarter-degree spatial resolution. Spatial averaging and quality assurance for the satellite data were the same as those in our previous study (Kondo et al 2015, Ichii et al 2017). A global map of forest classes for 2000 (Schulze et al 2019) was used for classifying forest/non-forest ecosystems. Rising [CO₂] compared with the level at year 2000 was determined based on a one-degree spatial resolution, where mean annual [CO₂] was calculated from CT2015. Assuming that fluxes based on the upscaling of the eddy covariance observations (Kondo et al 2015) already includes the CO₂ fertilization effect, finally, the CO₂ fertilization effect in the fluxes (${F}_{{{\rm{CO}}}_{2}};$ g C m⁻² yr⁻¹, ${\rm{E}}{{\rm{T}}}_{{{\rm{CO}}}_{2}};$ mm yr⁻¹) can be defined as follows:

$$\begin{eqnarray}&&\begin{array}{l}{F}_{{{\rm{CO}}}_{2}}=Fc* {f}_{{\rm{C}}3}-Fc* {f}_{{\rm{C}}3}* \\ \,\,\,\,\,\times \,(100.0/(100.0+{F}_{ \% {\rm{C}}}* {\rm{\Delta }}{{\rm{CO}}}_{2})),\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\rm{E}}{{\rm{T}}}_{{{\rm{CO}}}_{2}}=T\ast {f}_{{\rm{C}}3}-T\ast {f}_{{\rm{C}}3}\ast \\ \,\times \,(100.0/(100.0+{F}_{{\rm{ \% }}{\rm{E}}}\ast {\rm{\Delta }}{{\rm{CO}}}_{2})),\end{array}\end{eqnarray} \tag{ 2 }$$

where Fc and T are annual GPP (g C m⁻² yr⁻¹) and transpiration (mm yr⁻¹), respectively, in each year derived from empirical data-driven upscaling EC fluxes using the support vector regression (SVR) (Kondo et al 2015). The upscaled global fluxes based on SVR were reprocessed using latest MODIS data (collection 6). The first term of the right side is fluxes of C3 vegetation that includes the CO₂ fertilization effect, and the second term is fluxes of C3 vegetation that does not include the CO₂ fertilization effect. The denominator of the second term increases with increasing [CO₂] since 2000, and the effect of the CO₂ fertilization is subtracted from the second term according to the change of flux with respect to the increasing [CO₂] in each grid. Transpiration was partitioned from the upscaled ET based on an empirical method (Wei et al 2017). Satellite-derived LAI and land cover data (MOD12Q1) were used for calculating the fraction of transpiration from ET at an eight-day interval. f_C3 was the fraction of C3 vegetation within a grid based on the global map for C4 vegetation percentage (Still et al 2009). For multiplying C3 fraction, our estimates of the global CO₂ fertilization effect only consider the C3 vegetation. F_%C and F_%E were changes in Fc and T, respectively, per change in [CO₂] (% ppm⁻¹) derived from the random forest regressions tuned for the site-scale analysis (shown in above paragraph), and ΔCO₂ was the change in [CO₂] in each year since 2000 (ppm yr⁻¹). An ensemble mean of the 20 different random forest regressions was used to estimate the CO₂ fertilization effect and its standard deviation was used for the interval. The declining response of the CO₂ fertilization effect (figures 1(c), (f)) was considered in the upscaling, using a regression shown in figures 1(c), (f).

Zoom In Zoom Out Reset image size

For estimating potential uncertainties in the upscaled CO₂ fertilization effect associated with the input data, we compared the upscaled results by replacing the input data. First, we replaced the input global climate data from the CRU/NCEP data to ERA5 reanalysis data (DOI:10.24381/cds.68d2bb30). Second, we replaced the input MODIS-based LAI data to the Global Inventory Modeling and Mapping Studies (GIMMS) LAI3g product (Zhu et al 2013), which was based on the advanced very high resolution radiometer (AVHRR). The CRU data are based on station observations, which is less uncertain than the reanalysis data. The MODIS-based LAI is generally more accurate than LAI by the AVHRR. Since the upscaled CO₂ fertilization effects did not differ in terms of the different input data (table S4), we showed the upscaled results based on the combination of the CRU/NCEP and MODIS LAI data as our best guess. The small sensitivity to the input data could be because the current upscaling only used the annual aggregated variables (figure S6); thus, biases at the short timescales in the input data did not significantly propagated to the upscaled estimates.