Atmospheric CO 2 concentrations [CO 2 ] have risen at a rate of 2.16±0.09 ppm yr -1 in recent decades due to human activities (Le Quéré et al 2018) and will continue to increase unless emission reductions occur (Anderson et al 2019). Rising [CO 2 ] has increased photosynthesis at the leaf level (Norby et al 2005, Wenzel et al 2016) and gross primary productivity (GPP) at the ecosystem scale (Wenzel et al 2016). This process, known as the CO 2 fertilization effect, arises because the current [CO 2 ] is too low to saturate the carboxylation in the leaf, and thus limits photosynthesis. Therefore, the CO 2 fertilization effect is considered a negative feedback process that may reduce [CO 2 ] and mitigate global warming (Booth et al 2012) through the potential enhancement of the land CO 2 sink. In addition to the carbon cycle, the hydrological cycle has similarly been impacted by the increase in [CO 2 ] that reduces stomatal conductance (Keenan et al 2013, Frank et al 2015, Cheng et al 2017) and results in a decrease in continental evapotranspiration (ET) and an increase in water runoff (Betts et al 2007).
The magnitude of the CO 2 fertilization effect is not well understood because of the difficulty of direct measurements that can disentangle this effect from other processes (Norby et al 2005), limitation in satellite observations (de Kauwe et al 2016b), and the uncertainty in the ecosystem model (Friedlingstein et al 2014, Smith et al 2016). Free-air CO 2 enrichment (FACE) experiments can provide a direct estimate of the CO 2 fertilization effect by applying a step change in [CO 2 ] (Norby and Zak 2011) to the surrounding environment, but their application is expensive and laborious. Consequently, current experiments are located in easily accessible regions, such as young temperate forests and croplands (Long et al 2004, Körner 2009, Frank et al 2015). Furthermore, while FACE experiments have assessed the CO 2 fertilization effect under, for example, a doubling of [CO 2 ], current ecosystems are exposed to a gradually rising [CO 2 ] and their response may differ from that found in the FACE experiments. Although FACE data sets have improved vegetation models representing the ecosystem response to rising [CO 2 ] (Medlyn et al 2015), the mechanisms of the CO 2 fertilization effect incorporated into state-of-the-art earth system models has been a major source of uncertainty in climate projections through the internal feedbacks between photosynthesis and other related processes (e.g. respiration, biomass allocation, and mortality) (Booth et al 2012, Wenzel et al 2016, Friedlingstein et al 2014, Churkina et al 2009).
Here, we present a new quantification of the CO 2 fertilization effect at the ecosystem scale during the past two decades based on eddy-covariance (EC) flux measurements across arctic, boreal, temperate, and tropical ecosystems at the global scale. Although previously the CO 2 fertilization effect was inferentially evaluated for a smaller and shorter data set (Keenan et al 2013), we expanded the analysis by using a 104-flux tower site (770 site years) data set (Ichii et al 2017, Pastorello et al 2017), consisting of forests, grasslands, savanna, shrub, wetlands, tundra, and cropland ecosystem (figure S1 is available online at stacks.iop.org/ERL/15/084009/ mmedia, table S1). Our approach used the global network of the direct observations, and quantified the sensitivity of GPP and ET to rising [CO 2 ] with multiple constraints to avoid confounding effects and artifacts. We emphasize that our approach avoids a direct use of decadal trends in the fluxes for estimating the CO 2 fertilization effect, which is mostly too small to be detected with common statistical analyses using observational data directly (Baldocchi et al 2018).
To isolate the CO 2 fertilization effect from other confounding effects, we developed a technique for deriving model parameters that are sensitive to the fertilization effect, using observed diurnal, seasonal, and interannual variations of carbon and water fluxes as well as climate drivers and leaf area index (LAI) (figure S2). The method is similar to the so-called onepoint method commonly used in the leaf scale analyses (de Kauwe et al 2016a) that is based on the idea that light-saturated GPP is regulated by CO 2 concentration. The responses of water vapor flux and GPP to light, humidity, wind speed, and [CO 2 ] were constrained with direct observations during the past two decades (table S3) based on semi-continuous EC flux measurements across arctic, boreal, temperate, and tropical ecosystems at the global scale. The photosynthetic responses to light and [CO 2 ] were constrained by optimizing the biogeochemical model (A1-1 in supplemental material); stomatal responses to humidity and [CO 2 ] were constrained by the stomatal conductance model (A1-2 in supplemental material); and the responses to wind speed were constrained by the aerodynamic conductance model (A1-4 in supplemental material). The method has been previously applied for leaf-scale (Reed et al 1976), and extended to the canopy scale (Wang et al 2007, Ueyama et al 2016, 2018). After accounting for other confounding effects, such as radiative transfer (A1-3 in supplemental material), boundary layer conductance (A1-4 in supplemental material), and evaporation (A1-6 in supplemental material), we isolated responses of GPP and ET to current gradually rising [CO 2 ] by fitting the fluxes into the CO 2 demand and supply functions for photosynthesis in the model after considering other environmental effects (e.g. change in aerodynamic and canopy conductance).
We used a 104 flux tower sites (770 site years) (Ichii et 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 2 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 2 fertilization effect (Ueyama et al 2016(Ueyama et al , 2018)) 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 2 ], 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 2 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 2 ] measurements at the observation sites, and for consistency among the optimization and CO 2 fertilization experiment, S1). If time-variant LAI was not available, LAI was derived from MODIS MOD15A2H after smoothing the signal using TIME-SAT (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 2 ]. Parameters that are sensitive to CO 2 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 2 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 2 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 2 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 -2 s -1 and 34.6 W m -2 , respectively; slopes and R 2 between observations and the parameterized model also supported the validity of the optimization (figure S5).
Since the actual parameters for CO 2 fertilization processes (CO 2 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 2 ] for a given day. The CO 2 fertilization effect was calculated at the halfhourly 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 -2 s -1 . The quantified CO 2 fertilization effect represents the changes in GPP and gs by changes in [CO 2 ] whilst allowing environmental conditions (e.g. We separated two processes related to the changes in GPP and gs, namely the passive response by a hyperbolic response of photosynthesis to [CO 2 ] 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 2 fertilization experiment was conducted using median seasonality in the parameters instead of daily parameters in each year. CO 2 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 2 ], were available.
We tested how the EC data constrained the CO 2 fertilization effect using data from 14 sites selected from the FLUXNET2015 data set (table S2). As a nonconstrained 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 2 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 2 fertilization for GPP and gs.
For evaluating the CO 2 fertilization effect at the global scale, the effect examined at the site level was upscaled using a random forest regression, which is a machinelearning based regression. Changes in fluxes per change in [CO 2 ] (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 nonforest), 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 2 score was 0.61±0.05 for the effects for GPP and 0.78±0.05 for gs.
The CO 2 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 2 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 2 ] 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) ). An ensemble mean of the 20 different random forest regressions was used to estimate the CO 2 fertilization effect and its standard deviation was used for the interval. The declining response of the CO 2 fertilization effect (figures 1(c), (f)) was considered in the upscaling, using a regression shown in figures 1(c), (f).
For estimating potential uncertainties in the upscaled CO 2 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 2 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.
On average across all biomes, GPP increased concurrently with rising [CO 2 ] at the rate of 0.273%±0.014% yr -1 (figure 1(a)) or 0.138%±0.007% ppm -1 (figures 1(d), (g)) from 1990 to 2014; plus-minus sign denotes 95% confidence interval. The change was greater for nonforest than forest sites (p<0.01; figures 1(a), (b)), consistent with FACE experiments which showed a higher increase in GPP at ecosystems with lower leaf area (Norby and Zak 2011). This was probably because light more regulated photosynthesis of forests that had greater leaf area than non-forests, and thus greater contributions of shade-leaves weakened the sensitivity to [CO 2 ] in forests than non-forests. Stomatal conductance decreased at a rate of -0.143%±0.012% yr -1 (figure 1(b)) or -0.073%±0.006% ppm -1 (figure 1(e)) over the same time period. The decrease in gs was significantly greater in forests than non-forests (p<0.01; figures 1(d), (e)), because greater increase in GPP by rising [CO 2 ] mitigated decrease in gs in non-forests. Given the increase in GPP and the decrease in gs, intrinsic water use efficiency (iWUE; defined as GPP/gs) increased at the rate of 0.400%±0.011% yr -1 (figure 1(c)) or 0.202%±0.006% ppm -1 (figure 1(f)). The magnitude of the CO 2 fertilization effect (figures 1(d)-(f)) was roughly comparable to the values observed from previous FACE assessments in a temperate zone (0.13% ppm -1 for net ecosystem productivity (NEP); Norby et al 2005), and 0.2%-0.3% ppm -1 for iWUE; Frank et al 2015) and global-scale modeling that was constrained with [CO 2 ] growth rate for high-latitudes (0.13%±0.03% ppm -1 for NEP) and the extratropics (0.11% ±0.03% ppm -1 for NEP) (Wenzel et al 2016).
Importantly, the CO 2 fertilization effect for GPP and iWUE slightly declined with rising [CO 2 ] (p0.01) (figures 1(g), (i)). This decline of the CO 2 fertilization effect could be caused by two different processes: the passive response by a hyperbolic response of photosynthesis to [CO 2 ] and/or ecophysiological acclimation, known as down regulation. To isolate the two processes, we quantified the ecophysiological acclimation, assuming that the ecophysiological acclimation occurred with interannual/decadal changes in the ecophysiological parameters (Vc max25 , J max25 , m bb , and b bb ; shown in method). The simulations suggested that the down regulation mechanism played a marginal role on the estimated magnitudes in the CO 2 fertilization effect. This in turn indicates that the canopy-integrated parameters did not change significantly during the study period. The magnitude of the change from ecophysiological acclimation was one order of magnitude smaller (figure 2) than the passive response to rising [CO 2 ] (Katul et al 2000) with unchanged parameters (figures 1(g)-(i)). Nevertheless, the ecophysiological acclimation rather amplified the CO 2 fertilization effect in GPP at 54% of the ecosystems that we examined, and dampened the decrease in gs at 63% of the ecosystems (figure 2). This is possibly because a decrease in photosynthetic capacity of leaves was compensated by increasing LAI. Long-term field studies for partitioning change in leafscale parameters and LAI need to disentangle the compensation.
The spatial variability of the isolated CO 2 fertilization effect on GPP and gs was explained environmental drivers, such as mean annual air temperature, growing season length and temperature, and land cover type (figure S6; R 2 =0.59±0.05 for GPP and R 2 =0.72±0.11 for gs). The increase in GPP due to rising [CO 2 ] was greater in warmer regions than in colder regions (R 2 =0.29, p<0.01; figure S7). This was likely because rising [CO 2 ] could have effectively mitigated photorespiration under higher temperatures (Cernusak et al 2013). The fractional decrease in gs tended to be greater in ecosystems in colder rather than warmer climates (R 2 =0.14; p<0.01), probably because gs is related to GPP (Ball et al 1987), and thus the greater fertilization effect on GPP alleviated the decrease in gs in warmer ecosystems.
The upscaled CO 2 fertilization effect of the global C3 vegetation accounted for a substantial impact on the current increasing trend in global GPP (figure 3(a)). Based on the upscaled CO 2 fertilization effect using random forest regression (figure S6) and a data-driven global GPP estimation (Kondo et al 2015), we found that rising [CO 2 ] (29.0 ppm between 2000 and 2014; Peters et al 2007) enhanced GPP at a rate of 0.08% ppm -1 at the global scale, which was equivalent to an increase of 0.16% yr -1 , or 0.23 Pg C yr -2 . This enhancement is approximately 60% of the current increasing trend in global GPP (Kondo et al 2015) (0.38 Pg C yr -2 ; red line in figure 3(a)) and other studies (0.59±0.12 Pg C yr -2 ; Cheng et al 2017), indicating that the CO 2 fertilization is the most important process driving the current increase in the global GPP. Enhanced GPP from CO 2 fertilization after the baseline year 2000 is, on average, 1.2% of global GPP (152.8 Pg C yr -1 ), 12.4 g C m -2 yr -1 or 1.8 Pg C yr -1 at During the study period, the global GPP increased (p<0.01), but diminished the trend if the CO 2 fertilization effect estimated in this study was excluded (figure 3(a)). Note that our estimates of the global CO 2 fertilization effect were only quantified for C3 vegetation, neglecting contributions by C4 vegetation, and thus the estimated magnitude was likely underestimated.
Our analysis showed that rising [CO 2 ] dampened an increase in the global ET (figure 3(b)). The global ET was estimated to increase at a rate of 1.04 mm yr -2 . This increase resulted from a positive balance between the negative effect of decreased gs due to rising Growing number and further long-term monitoring of eddy covariance observations will further reduce the uncertainties in the CO 2 fertilization effect. Limited data in the current analysis is available for data whose temporal extent longer than nine years and for tropics where the CO 2 fertilization effect estimated to be large, resulting in biased estimates in the global upscaling. This is the major shortcoming in our estimates, but will be minimized with future FLUXNET activities by further covering wide range of [CO 2 ], climate, and ecosystem types.
The CO 2 fertilization effect, the increase in GPP and decrease in gs, should be occurring in global ecosystems at present, ranging from the tropics to the arctic. The
Res. Lett. 15 (2020) 084009