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

Rising atmospheric CO2 concentration ([CO2]) enhances photosynthesis and reduces transpiration at the leaf, ecosystem, and global scale via the CO2 fertilization effect. The CO2 fertilization effect is among the most important processes for predicting the terrestrial carbon budget and future climate, yet it has been elusive to quantify. For evaluating the CO2 fertilization effect on land photosynthesis and transpiration, we developed a technique that isolated this effect from other confounding effects, such as changes in climate, using a noisy time series of observed land-atmosphere CO2 and water vapor exchange. Here, we evaluate the magnitude of this effect from 2000 to 2014 globally based on constraint optimization of gross primary productivity (GPP) and evapotranspiration in a canopy photosynthesis model over 104 global eddy-covariance stations. We found a consistent increase of GPP (0.138 ± 0.007% ppm−1; percentile per rising ppm of [CO2]) and a concomitant decrease in transpiration (−0.073% ± 0.006% ppm−1) due to rising [CO2]. Enhanced GPP from CO2 fertilization after the baseline year 2000 is, on average, 1.2% of global GPP, 12.4 g C m−2 yr−1 or 1.8 Pg C yr−1 at the years from 2001 to 2014. Our result demonstrates that the current increase in [CO2] could potentially explain the recent land CO2 sink at the global scale.


Introduction
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  (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 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 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).

Method
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).

Data
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 2 fertilization effect is expected in C4 ecosystems than C3-dominated ecosystems (Collatz et al 1992).

Model and optimization
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. The model contains coupled effects among the following sub- 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 2 concentration, whereas RuBP-limited photosynthesis is sensitive to light. Stomatal conductance is modeled based on a semi-empirical relationship to photosynthetic rate, CO 2 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 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)  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. climate and LAI) to change over time (Ueyama et al 2016). Here, we used [CO 2 ] in the year 2000 as a baseline (369.55 ppm of the annual CO 2 concentration at Mauna Loa) (Le Quéré et al 2018) and evaluate the effects of changing [CO 2 ] 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 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.

Upscaling
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) 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 2 ] compared with the level at year 2000 was determined based on a one-degree spatial resolution, where mean annual [CO 2 ] was calculated from CT2015. Assuming that fluxes based on the upscaling of the eddy covariance observations (Kondo et al 2015) already includes the CO 2 fertilization effect, finally, the CO 2 fertilization effect in the fluxes (F ; CO 2 g C m −2 yr −1 , ET ; CO 2 mm yr −1 ) can be defined as follows: where Fc and T are annual GPP (g C m −2 yr −1 ) and transpiration (mm yr −1 ), 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 2 fertilization effect, and the second term is fluxes of C3 vegetation that does not include the CO 2 fertilization effect. in each year since 2000 (ppm yr −1 ). 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.

Results and discussion
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 Figure 1. Inferred fractional changes in gross primary productivity (GPP) (a), canopy-integrated stomatal conductance (gs) (b), and intrinsic water use efficiency (iWUE; GPP/gs) (c) associated with rising atmospheric CO 2 concentration ([CO 2 ]). The CO 2 fertilization effect was quantified with changes in [CO 2 ] since the year 2000, which was taken as baseline (ΔCO 2 ). The fractional changes were shown as the ratio of GPP, gs, or iWUE calculated with and without the CO 2 fertilization effect. Analysis was performed based on the coupled photosynthesis and stomatal conductance model, with optimized parameters for each day and each site. Dots and bars represent the mean and range of uncertainties with different formulations in the photosynthesis models, respectively. Relationship between fractional change in GPP (d), gs (e), and iWUE (f) associated with rising [CO 2 ] and change in annual mean [CO 2 ]. Relationship between fractional change in GPP (g), gs (h), and iWUE (i) per unit of [CO 2 ] change and change in annual mean [CO 2 ]. Square dots in (g)-(i) represent averaged values for ΔCO 2 bins of two ppm. The statistically significant regressions (p0.01) are shown in (d)-(g), (i).
Importantly, the CO 2 fertilization effect for GPP and iWUE slightly declined with rising [CO 2 ] (p0.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) 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 Figure 2. Probability density distributions of the difference between the CO 2 fertilization effects in GPP (a), canopy-integrated stomatal conductance (gs) (b), and intrinsic water use efficiency (iWUE) (c), estimated using daily optimized parameters and using median seasonal variation of the parameters. The analysis was done for 31 eddy covariance sites having more than four years of data. Estimates using the median seasonality remove year-to-year and decadal shift in physiological adjustment; thus, the difference represents change in CO 2 fertilization effect associated with the physiological adjustment. A positive value represents an increase in GPP and iWUE or dampened decrease in gs associated with physiological adjustment. the years from 2001 to 2014. This increase was surprisingly large, because the current global land sink was estimated to be 3.5±1.0 Pg C yr −1 during 2008-2017 (Le Quéré et al 2018). The tropical forests benefited substantially from the CO 2 fertilization effect (figure 3(e)), due to greater effects per rising [CO 2 ] (figure 3(c)) and greater GPP in tropical forests. 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 [CO 2 ] and the positive effect of the enhanced atmospheric evaporative demand induced by global warming (Katul andNovick 2009, Cheng et al 2017). In practice, we inferred that without rising [CO 2 ], the global ET would have increased at a rate of 1.24 mm yr −2 because the regulation of gs by rising [CO 2 ] dampened global transpiration (Wullschleger et al 2002), resulting in the reduction of ET by a rate of 0.20 mm yr −2 ( figure 3(b)). for ET, whereas those for absolute change (unit for g C m −2 yr −1 for GPP and mm yr −1 for ET) are shown in (e) for GPP and (f) for ET. The CO 2 fertilization effect was quantified, with the atmospheric CO 2 concentration in 2000 as a baseline. Note that upscaled ET was the sum of evaporation and transpiration, whereas the CO 2 fertilization effect occurs only transpiration; the magnitude of the change in ET (b), (d), (f) was estimated after partitioning transpiration from ET. The shadows represent the range of uncertainties associated with different model parameterization, energy imbalance, intercepted evaporation, and upscaling.
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.

Conclusion
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 quantified effect was comparable to results from previous field experiments in temperate regions (Norby et al 2005, Frank et al 2015, reinforcing the validity of our approach. Our study enabled spatial upscaling from available flux tower data to provide a global understanding of tropical, boreal, and arctic regions, where direct measurements are scarce or simply not available, and gives insight into the possible mechanisms of the current land CO 2 sink and decadal trends in ET. By constraining a simple model with observed big data, the direct CO 2 fertilization effect was estimated globally in this study, but other indirect effects, such as, changes in LAI and extra-growth by relaxed drought owing to the

Acknowledgments
This study was supported by the JSPS KAKENHI, Grant Number 16K12588. Databases of eddy covariance measurements (AsiaFlux, KoFlux, European Fluxes Database Cluster, FFPRI, NIAES, and the NSF Arctic Observatory Network) facilitated this study. This work used eddy covariance data, FLUXNET2015, acquired and shared by the FLUXNET community, including those networks: AmeriFlux, AsiaFlux, NECC, OzFlux TERN, TCOS-Siberia, and USCCC. The ERA-Interim reanalysis data provided by ECMWF and processed by LSCE. The FLUXNET eddy covariance data processing and harmonization was carried out by the European Fluxes Database Cluster, AmeriFlux Management Project, and Fluxdata project of FLUXNET, with the support of CDIAC and ICOS Ecosystem Thematic Center, and the OzFlux, ChinaFlux and AsiaFlux offices. LBM acknowledges the support of the 'RUDN University program 5-100'. The GIMMS LAI3g product was provided by Professor R B Myneni.

Data availability
All flux data are available from respective databases or by authors upon request. (https://doi.org/10.25412/iop. 11559480.v1) The program codes of the iBLM-EC are available from the author's web-site (http://atmenv. envi.osakafu-u.ac.jp/staff/ueyama/softwares/) or a supplemental material from the IOP Publishing Figshare repository. The inferred daily ecophysiological parameters, annual CO 2 fertilization effects, and the gridded data for the global CO 2 fertilization effects are available through a supplemental material from the IOP Publishing Figshare repository.