Global decadal variability of plant carbon isotope discrimination and its link to gross primary production.

Carbon isotope discrimination (Δ13 C) in C3 woody plants is a key variable for the study of photosynthesis. Yet how Δ13 C varies at decadal scales, and across regions, and how it is related to gross primary production (GPP), are still incompletely understood. Here we address these questions by implementing a new Δ13 C modelling capability in the land-surface model JULES incorporating both photorespiratory and mesophyll-conductance fractionations. We test the ability of four leaf-internal CO2 concentration models embedded in JULES to reproduce leaf and tree-ring (TR) carbon isotopic data. We show that all the tested models tend to overestimate average Δ13 C values, and to underestimate interannual variability in Δ13 C. This is likely because they ignore the effects of soil water stress on stomatal behavior. Variations in post-photosynthetic isotopic fractionations across species, sites and years, may also partly explain the discrepancies between predicted and TR-derived Δ13 C values. Nonetheless, the "least-cost" (Prentice) model shows the lowest biases with the isotopic measurements, and lead to improved predictions of canopy-level carbon and water fluxes. Overall, modelled Δ13 C trends vary strongly between regions during the recent (1979-2016) historical period but stay nearly constant when averaged over the globe. Photorespiratory and mesophyll effects modulate the simulated global Δ13 C trend by 0.0015 ± 0.005‰ and -0.0006 ± 0.001‰ ppm-1 , respectively. These predictions contrast with previous findings based on atmospheric carbon isotope measurements. Predicted Δ13 C and GPP tend to be negatively correlated in wet-humid and cold regions, and in tropical African forests, but positively related elsewhere. The negative correlation between Δ13 C and GPP is partly due to the strong dominant influences of temperature on GPP and vapor pressure deficit on Δ13 C in those forests. Our results demonstrate that the combined analysis of Δ13 C and GPP can help understand the drivers of photosynthesis changes in different climatic regions.

INTRODUCTION

Plant tissues have a 13C‐depleted isotopic signature relative to atmospheric CO2 because of the slower diffusion of heavier 13CO2 molecules through the stomata, and the preferential fixation of 12CO2 by the enzyme RuBisCO during carboxylation. This isotopic fractionation is well‐known as the discrimination against 13C, and denoted Δ13C (Park & Epstein, 1960). Despite decades of research, however, there is still incomplete knowledge of how and why Δ13C in C3 woody plants varies over decadal timescales, and across regions.

Studies investigating short‐term (seasonal to interannual) variations in Δ13C and their spatial patterns using stable carbon isotope ratios (δ13C) measured in C3 leaves have shown that changes in Δ13C depend on environmental conditions, including soil moisture and annual precipitation (Diefendorf et al., 2010; Kohn, 2010), atmospheric vapor pressure deficit (Lloyd & Farquhar, 1994; Wang, Prentice, Keenan, et al., 2017), temperature, atmospheric pressure (Cornwell et al., 2018), and thus elevation (Korner et al., 1988, 1991; Wang, Prentice, Davis, et al., 2017)—see Cernusak et al. (2013) for a review. Other studies using δ13C measured in tree rings (TRs) have suggested that Δ13C in C3 woody plants increases with rising atmospheric CO2 (Ehleringer & Cerling, 1995; Feng & Epstein, 1995; McCarroll et al., 2009; Treydte et al., 2001). Experimental studies investigating changes in Δ13C of (well‐watered) herbaceous C3 angiosperm plants growing at different CO2 levels (Schubert & Jahren, 2012, 2015) have also shown this pattern. In principle, Δ13C might increase with rising CO2 in mesic or well‐watered sites when stomata are fully open and thus the level of CO2 inside the leaf (c i) is high, reflecting the RuBisCO enzyme's preference for 12CO2 (Farquhar et al., 1982). However, even with rising CO2, Δ13C might decrease in dry regions when stomatal conductance is reduced and c i is low (Farquhar et al., 1982).

The actual sensitivity of Δ13C to changes in CO2 continues to be debated. For example Sheldon et al. (2020) found no significant change in Δ13C values derived from a large global data set of leaf measurements from woody gymnosperms for the post‐1850 period. They suggested that fundamental differences in Δ13C responses to CO2 between plant functional groups might explain the apparent discrepancies between experimental and historical studies. Consistent with these findings, Hare and Lavergne (2021) used a model‐data fusion approach to show that the Δ13C in woody gymnosperms is less sensitive to CO2 than that in angiosperms. It has also been suggested that decadal‐scale plant physiological responses, including adjustments of stomatal size and density, could counteract the positive effect of CO2 on Δ13C values in woody C3 plants (Rayback et al., 2020; Saurer et al., 2004; Stein et al., 2020).

One study exploring 13CO2 budgets using atmospheric measurements and a box model (Keeling et al., 2017) estimated that plant Δ13C should have increased globally by 0.014 ± 0.007‰ ppm–1 over the period 1978–2014, in order to explain the magnitude of the decreasing trend observed in atmospheric δ13CO2 (attenuation of the Suess effect: Keeling, 1979). Suggested reasons for such an increase were the impact of CO2 on photorespiratory fractionation and (more importantly) mesophyll‐conductance fractionation. Their results imply that despite differing temporal variations in Δ13C among sites, Δ13C should have shown a global increase, when integrated over all ecosystems. However, the global Δ13C trend suggested by this study has not been independently validated.

In addition to the global Δ13C trend, the degree to which Δ13C can be used as proxy for interannual variations in gross primary production (GPP) is still unclear. During large‐scale drought, when continental‐scale net carbon uptake is low, so too is Δ13C, leading to a positive relationship between Δ13C and GPP (Peters et al., 2018; Randerson et al., 2002). This relationship is consistent with Belmecheri et al. (2014) findings in the (mesic) Harvard Forest using TR and eddy‐covariance (EC) flux measurements. However, other studies not only reported positive relationships between Δ13C and TR growth in relatively dry regions, but also negative ones in wet and/or cold regions (del Castillo et al., 2014; Shestakova et al., 2019; Voelker et al., 2014). Because TR growth is an indicator of GPP (Babst et al., 2014), their results suggest that even the sign of the GPP‐Δ13C relationship could vary with environment.

Here we explore mean values and trends in global Δ13C of C3 woody plants over the 1979–2016 period using the Joint UK Land Environment Simulator (JULES) model (version 5.6), newly equipped with a carbon isotopic modelling capability. We first evaluate the model at the local scale, using a large network of δ13C data from leaves and TRs to assess which c i model in JULES shows the best agreement with observations. We then analyze the predicted Δ13C values and trends from the global historical simulations and quantify the relative contributions to the global trend. We also investigate the sign, magnitude, and drivers of the relationship between Δ13C and GPP across regions. We aim to answer three questions: (1) Is the global increase in plant Δ13C suggested by Keeling et al. (2017) reproduced by JULES? (2) What are the contributions from photorespiratory and mesophyll effects to the global Δ13C trend? and (3) Is the correlation between Δ13C and GPP driven by environmental conditions, as suggested by TR studies?

MATERIALS AND METHODS

Description of the JULES model

JULES is the land surface component of the UK Earth System Model (Sellar et al., 2019) simulating simultaneously the fluxes of carbon, water and energy between the land surface and the atmosphere (https://jules.jchmr.org). Vegetation in the model is represented as plant functional types (PFTs), which differ in their biochemical and biophysical properties. A full description of the model is available in previous publications (Best et al., 2011; Clark et al., 2011; Harper et al., 2016).

Leaf‐level photosynthesis and soil water stress factor

The potential assimilation rate (A p, μmol m−2 s−1) for C3 plants is predicted using the standard Farquhar et al. (1980) model as the minimum of two limiting factors (electron transport A J and carboxylation rate A C, see Text S1) minus the dark mitochondrial respiration (R d; μmol m−2 s−1). A C is determined by the maximum carboxylation rate (V cmax, μmol m−2 s−1), which regulates the RuBisCO enzymatic capacity for carbon fixation, and A J depends on light (with a maximum value determined by the maximum electron‐transport rate, J max, μmol m−2 s−1). Both limiting factors also depend on the leaf intercellular partial pressure of CO2 (c i, Pa) and on the CO2 compensation point in the absence of R d (Γ*, Pa). R d is proportional to V cmax but varies among PFTs through the parameter f d (Harper et al., 2016). The net assimilation rate (A n) resulting from biochemical limitations of photosynthesis due to soil water stress is then estimated as: W is the smoothed minimum of A C and A J. β soil is a dimensionless soil water stress factor (β soil) ranging between 1 (well‐watered plants) and 0 (no available water) and determined based on the soil characteristics (texture and type) following the Van Genuchten (1980) model.

β (unitless) is calculated in each grid box for each soil layer k as the function of volumetric soil water content in each layer of the root zone (, m3 m−3) using: with

is the critical above which plants are unaffected by water stress (defined as the field capacity, m3 m−3) and is the wilting below which water stress is at its maximum (or permanent wilting point, m3 m−3) in each soil layer k. is the water content at which plants start to become water stressed in each soil layer k. p 0 is a parameter governing the threshold at which the plant starts to experience water stress due to lack of water in the soil (ranging between 0 and 1). A p 0 value equal to 0.4 is used here to delay the onset of soil water stress, consistently with recent studies (Harper et al., 2021; Williams et al., 2019). The overall β soil is calculated based on the water stress in each layer estimated from Equation ((2a), (2b)) and the fraction of root mass in each soil layer k (r k): where n soil is the number of soil layers (equal to 14 in this study).

Leaf intercellular CO2 models tested in JULES

The leaf‐level stomatal conductance for carbon (g sc; mol m−2 s−1) is related to A n from Equation (1) via the CO2 diffusion equation as: where c a and c i are the ambient and leaf‐intercellular CO2 (in μmol mol−1 here), respectively.

Two representations of c i have been incorporated into JULES based on the Jacobs (1994) and Medlyn et al. (2011) stomatal conductance models. The Jacobs (1994) model originally implemented into JULES (Best et al., 2011) relates c i/c a to leaf humidity deficit as: d q and are the specific and critical humidity deficit at the leaf surface (kg kg−1). f 0 is the PFT‐dependent c i/c a at the leaf‐specific humidity deficit.

The Medlyn et al. (2011) model incorporated and validated in JULES version 5.5 (Oliver et al., 2018) assumed that stomatal aperture is regulated to maximize carbon gain while simultaneously minimizing water loss so that: D (kPa) is the leaf‐to‐air vapor pressure deficit and g 1 (kPa0.5) is a PFT‐dependent fitted parameter representing the sensitivity of g sc to A, that is, the normalized water use efficiency.

Two other stomatal models have been incorporated into JULES version 5.6 for this study, the Leuning (1995) and the Prentice et al. (2014) models. The Leuning (1995) model is a modified version of the simple Ball et al. (1987) model and is widely used in land surface modelling: g 0 (μmol m−2 s−1), D 0 (kPa), a 1 (unitless) are empirically fitted parameters representing the residual stomatal conductance and the sensitivity of stomata to changes in D and A, respectively.

The Prentice et al. (2014) model is based on the least‐cost optimality hypothesis which assumes that leaves minimize the summed unit costs of transpiration and carboxylation so that: where β (unitless) is the ratio of cost factors for carboxylation and transpiration at 25°C, which may vary with changes in plant‐available soil water (Lavergne et al., 2020), but is set constant here because no mechanistic formulation for the soil water stress has been proposed yet. K is the effective Michaelis constant for Rubisco‐limited photosynthesis at a given partial pressure of O2 (Pa). η* (unitless) is the dynamic viscosity of water relative to its value at 25°C, assumed constant here (equal to 1). Both Γ* and K are calculated from their respective values at 25°C and 99.1 kPa (and activation energy) following Bernacchi et al. (2001).

All these models (Equations (5), (6), (7), (8a), (8b), (9a), (9b), (10), (11), (12), (13), (14), (15a), (15b), (15c), (16), (17)) implicitly assume infinite mesophyll conductance (g m) and therefore that the ratio of chloroplastic (c c) to ambient (c a) CO2 (c c/c a) is equal to c i/c a. This is an oversimplification as g m is low enough to cause a drawdown from c i to c c (Flexas et al., 2012), leading to lower c c/c a compared to c i/c a. An alternative formulation for the Prentice model assuming finite g m has been proposed by Wang, Prentice, Keenan, et al. (2017) and is also tested here: where β (unitless) is equivalent to β c assuming finite g m. and K c are equivalent to Γ* and K assuming finite g m calculated following Bernacchi et al. (2002). The g sc/g m ratio is assumed constant here.

Upscaled canopy fluxes of carbon and water

The leaf‐level A n, R d and g sc in each layer of the canopy (i) are scaled by the fraction of leaf area index (LAI, m2 m−2) for sunlit and shaded leaves and summed over all model layers (n, equal to 10 here) to obtain canopy‐scale estimates A can, R dc and G sc as: where X is the variable considered and, f sun,i and f shade,i are the fractions of sunlit and shaded leaves for each canopy layer i, respectively.

The GPP (gC m−2 s−1) is then estimated as:

The factor 0.012 converts from mol CO2 m−2 s−1 to kg C m−2 s−1.

The surface latent heat flux (LE, W m−2) is calculated using the humidity gradient between the surface and atmospheric reference height z 1 following Cox et al. (1999): is the factor determined from the proportions of canopy evaporation, bare soil evaporation, transpiration by vegetation through the stomata, and sublimation from snow. L is the latent heat of vaporization of water (J kg−1), is the surface air density (kg m−3), r a is the aerodynamic resistance (s m−1), q sat (T s) is the saturated specific humidity at surface temperature T s (kg kg−1), and q 1 is the specific humidity at reference height z 1 (kg kg−1). LE is converted from energy (W m−2) into mass (evapotranspiration, ET in kgH2O m−2 s−1) by dividing LE with L.

The transpiration flux (T r, kgH2O m−2 s−1) is calculated for the vegetated fraction of the grid box f v (unitless) as: f a is the wet fraction of the canopy (unitless) and r c is the canopy resistance (s m−1) inversely related to the canopy g sw (G sw, m s−1). G sw is related to G sc following Equation (4) as: with R the ideal gas constant (J K−1 mol−1).

Implementation of stable carbon isotopes into JULES

Assuming infinite boundary layer conductance and negligible fractionation during mitochondrial respiration (Evans & von Caemmerer, 2013), Δ13C in C3 plants depends on the isotope fractionations due to CO2 diffusion across the stomata pore (a = 4.4‰) and the mesophyll cell (a m = 1.8‰), and due to effective RuBisCO carboxylation (b = 30‰) and photorespiration (f = 12 ± 4‰) following Farquhar et al. (1982):

Assuming infinite g m, Equation (15a) can be simplified as: where  = 28‰ following Ubierna and Farquhar (2014) to account for the missing mesophyll term. Finally, assuming no fractionation during photorespiration, Equation (15b) is simplified: where b′ = 27‰ following Farquhar et al. (1982) to account for the missing mesophyll and photorespiration terms. We incorporated these three variants of the discrimination model into JULES version 5.6. Note that for Equation (15a), only the Prentice stomatal model (Equation (9a), (9b)) was used because the other three models do not predict c c.

Model configurations, simulations, and evaluation with observations

Model setup: Configurations and environmental drivers

We ran the model on the NERC JASMIN platform (http://www.jasmin.ac.uk/) using the Rose/Cylc suite control system. We used the suite control setup to initialize, configure, spin up, and run the historical simulations. We used a JULES configuration (Rose suite u‐bx886) with 13 surface tiles consisting of nine PFTs (tropical broadleaf evergreen trees, temperate broadleaf evergreen trees, broadleaf deciduous trees, needle‐leaf evergreen trees, needle‐leaf deciduous trees, C3 grasses, C4 grasses, evergreen shrubs, and deciduous shrubs) following Harper et al. (2016) and four non‐vegetated surface types (urban, inland water, bare soil, and ice). The PFT‐dependent parameters values of the photosynthesis and c i models tested here for the five forest PFTs are reported in Table S1.

The model was forced by the WFDEI meteorological reanalysis data set at 0.5° × 0.5° spatial and 3 h temporal resolution over 1979–2016 (Weedon et al., 2014). All forcing variables were interpolated down to half‐hourly resolution. We ran JULES with prescribed annual mean atmospheric CO2 concentrations from NOAA ESRL (www.esrl.noaa.gov/gmd/ccgg/trends/). We used a fixed land cover mask based on the European Space Agency's Land Cover Climate Change Initiative global vegetation distribution (Poulter et al., 2015) and prescribed the soil properties using ancillary data based on the Van Genuchten (1980) model (Slevin, 2017). The model was spun up for a maximum of three cycles of 38 years length. The spin‐up simulations were performed starting with existing climatology provided by the meteorological data as initial condition, by iterative repeating the years of simulation (1979–2016) until the solutions converge on two in‐soil prognostic variables following these criteria: (1) for soil moisture: to within 1% of content, (2) for soil temperature: to within 0.1°K.

Model simulations at the point‐scale and evaluation with observations

We first evaluated the JULES model regarding its representation of the controls by leaf stomata of the coupled carbon and water cycles to determine which c i model should be used for the global simulations. We tested different combinations of the discrimination and stomatal models described above (see Table S2) running the model at the point‐scale where isotopic data for C3 woody plants were available and compared the resulting Δ13C predictions with the observational data set described below.

Simulations at sites with leaf or TR isotopic measurements

We used the δ13C data set used in Lavergne, Sandoval, et al. (2020) (referred to as ‘global network’, Figure S1a), which comprises leaf δ13C data from Diefendorf et al. (2010), Cornwell et al. (2018) and Sheldon et al. (2020) and TR δ13C data from Lavergne, Voelker, et al. (2020). We additionally used TR δ13C data derived from other recent compilations (Adams et al., 2020; Belmecheri et al., 2021; Mathias & Thomas, 2021). From an original compilation of 545 TR chronologies spanning the 1901–2016 period, we only selected those with at least 25 years of measurements over the 1979–2012 period (total of 151 chronologies). We reduced the analysis period to 1979–2012 (not 1979–2016) because only a few TR chronologies were available after 2012.

We then calculated ∆13C as: where δ13CO2 and δ13CPM (‰) are the stable isotopic compositions of the atmospheric CO2 and of the plant material considered (bulk leaf or α‐cellulose in TR), respectively. Yearly δ13CO2 values were derived from Graven et al. (2017). f post quantifies the sum of fractionations beyond those associated with the production of the primary photosynthetic assimilates—δ13C being higher in wood than in leaves (Cernusak et al., 2009). Because the leaf and TR data from the global network were not available at the same sites and, thus, the post‐photosynthetic fractionations could not be estimated for each TR site, we used a standard estimate (f post = 2.1 ± 1.2‰ between bulk leaf and α‐cellulose in TR; Belmecheri and Lavergne, 2020) assuming that the resulting ∆13C was equivalent to bulk leaf‐derived ∆13C. Thus, we considered that the fractionation effects beyond those occurring within the leaf were constant for all tree species and sites, and across years.

We then compared the leaf and TR‐derived Δ13C data with the predicted Δ13C values from the different combinations of the discrimination and stomatal models. We averaged the predicted Δ13C values over summer months (i.e., June–August for Northern Hemisphere and December–February for Southern Hemisphere). Alternatively, we tested weighting the predicted Δ13C for the PFT considered using the predicted GPP values for each year as:

For comparison of predicted values with the leaf‐derived Δ13C data set, we used Taylor diagrams (Taylor, 2001) to provide a statistical summary of the predictive skills of the model for the different simulations to reproduce the observational data set in terms of their Pearson's correlation coefficient (r), normalized standard deviation (SD) and centered root‐mean‐square error (RMSE). For comparison with the TR‐derived Δ13C timeseries, we estimated the Theil‐Sen's slope and interannual variability (IAV) of each model prediction and TR timeseries and compared the resulting values between each other at each site. The IAV was calculated as the SD of the difference between the Δ13C values and the Theil‐Sen's slope (i.e., the median of the slopes of all lines through pair of points; Sen, 1968).

Simulations at EC flux stations

We ran the model at EC flux stations where both carbon and water flux measurements and δ13C data in different parts of the plants were measured simultaneously for several years (referred to as ‘EC flux stations’, Figure S1b). The simulations at these sites allowed us to test the impact of the stomatal model selected on the predicted canopy fluxes. Using leaf and TR δ13C data at these sites we can also quantify the post‐photosynthetic fractionations and compare them to the standard value we used for the other TR sites.

We used the isotopic data set from Hemming et al. (2005) providing δ13C measurements in bulk leaf and stem at nine EC stations from the EU CarboEuropeFlux project in 2001 and/or 2002. We also exploited the published δ13C data set in bulk leaves for a few individual years (Guerrieri et al., 2016) and α‐cellulose in TR for the 1982–2012 period (Guerrieri et al., 2019) at eight AmeriFlux stations (see Table S3). We then estimated the apparent f post for each site and each year as the difference between the leaf and wood δ13C values. For the Hemming et al. (2005) data set, we corrected the obtained f post values by 1.1‰ to account for differences between bulk and α‐cellulose TR δ13C following Guerrieri et al. (2017) estimates.

Half‐hourly or hourly GPP and LE data were extracted from the FLUXNET‐2015 Tiers 1 product (GPP_NT_VUF_REF and LE_F_MDS, respectively) when available for sites studied for the CarboEuropeFlux (BE‐Bra, DE‐Hai, DE‐Tha, DK‐Sor FI‐Hyy, FR‐LBr, IT‐Col, It‐Lav, and NL‐Loo) and AmeriFlux (US‐MMS and US‐Ha1) stations (Table S3). We screened and filtered the data as follows: (1) we used data flagged as good, (2) we removed data errors (D < 0.01 kPa), (3) we removed nighttime data using the nighttime/daytime flag based on potential radiation when provided or using the photosynthetic photon flux density (PPFD) variable (PPFD >0 during daytime), and (4) we excluded non transpirational water fluxes by removing days with precipitation and the subsequent 24 h. We converted GPP into gC m−2 day−1 and estimated ET (in kg H2O m−2 day−1) from LE and air temperature. We then aggregated GPP and ET into means over the growing season (May–September) and compared the resulting observations with the simulations for the same period using Taylor diagrams.

Global simulations over 1979–2016

We ran the model at the global scale over the 1979–2016 period using the best stomatal model found with the point‐scale simulations (referred to as ‘global run’) and analyzed global mean values and trends in Δ13C predicted by JULES for the different versions of the discrimination model. For this global run, we assumed an isotopic fractionation due to Rubisco carboxylation of 30‰ for all discrimination models (Equations (15a), (15b), (15c)) to allow a fair intercomparison of the simulations. The predicted Δ13C for the five forest PFTs were weighted for each PFT using the GPP values and then averaged as in Equation (17) to produce one weighted‐Δ13C value per year and per grid‐point. In the following we use the term Δ13C when referring to weighted‐Δ13C. To explore the environmental drivers of Δ13C trends across the globe, we grouped the grid‐points for ranges of 5°C in annual air temperature (T air), 0.4 kPa in vapor pressure deficit (D) and 0.1 in soil water availability (β soil) and calculated the average of the Δ13C trends within each group (referred to ‘trend score’). We used a similar approach for exploring patterns of GPP trends across the globe.

To quantify the relative contributions of photorespiratory and mesophyll effects to the global Δ13C trend, we compared the simulation runs for the different discrimination models. To determine the sign and magnitude of the Δ13C‐GPP relationship across the globe and their links to environmental conditions, we calculated the Pearson's r between simulated annually averaged Δ13C and GPP at each grid point. We then grouped the resulting Δ13C‐GPP grid‐point correlations for the different ranges of annual T air, D and β soil values and calculated the average of the Pearson's r within each group (referred to ‘correlation score’). Finally, to investigate the relative environmental dependencies of Δ13C and GPP within each correlation score group, we performed multiple regression models of Δ13C and GPP against their common drivers (i.e., T air, D and atmospheric CO2) for the different groups.

RESULTS

Selection of the best configuration of JULES using a data‐model approach

Predicted mean Δ13C values: Comparisons with leaf isotopic measurements

The observed leaf‐derived Δ13C averaged 20.1 ± 2.3‰ (mean ± SD; n = 1558; Figure 1a) and was higher in tropical areas, northern Europe and the east coast of northern America, but lower in mountainous regions, mainland Australia and the west coast of northern America (Figure 1c). All Δ13C predictions derived from the four stomatal models tested in JULES tended to overestimate median Δ13C and underestimate the SD in Δ13C as observed in the leaf isotopic data set (Figure 1a,b). Using the simple version of the discrimination model, the Δ13C values were overestimated by 2.8 ± 0.8‰ but incorporating the photorespiratory effect reduced the difference by 0.4‰ (2.4 ± 0.7‰; Figure 1a). The incorporation of the mesophyll effect in the Prentice model resulted in a further reduction of the bias, leading to a difference of around 1.1‰ between observed and predicted Δ13C median values. Overall, using the Prentice model combined with the discrimination model incorporating both photorespiratory and mesophyll effects produced the highest Pearson's r and the lowest RMSE (Table 1; Figure 1b).

FIGURE 1

Impacts of stomatal and discrimination representations in JULES on the predicted Δ13C values as compared to the global network of leaf Δ13C measurements. Boxplots (a) and Taylor diagram (b). Leaf‐derived Δ13C values across the globe (c)

TABLE 1

Pearson's correlation coefficients (r) between leaf‐derived and predicted Δ13C values for each stomatal and discrimination model. All correlations are significant at p < .001. Summer months (i.e., June–August for Northern Hemisphere and December–February for Southern Hemisphere)

Stomatal model	Discrimination model	Summer months	Weighted by GPP
Jacobs	Simple	0.51	0.36
	Photorespiration	0.46	0.29
Leuning	Simple	0.31	0.36
	Photorespiration	0.29	0.27
Medlyn	Simple	0.51	0.28
	Photorespiration	0.49	0.24
Prentice	Simple	0.52	0.43
	Photorespiration	0.54	0.42
	Photorespiration + finite mesophyll	0.54	0.44

Predicted Δ13C trends and IAV: Comparisons with tree‐ring isotopic timeseries

Overall, TR Δ13C stayed relatively constant when considering all sites together (p = .097), but we found different patterns of variations when looking at the individual TR timeseries. While Δ13C in 95 out of 151 sites (63%) stayed nearly constant over 1979–2012, 29 sites showed a significant increasing trend (average 0.034 ± 0.006‰ year−1) and 27 depicted a negative trend (−0.037 ± 0.07‰ year−1). The sites with positive trends were mainly located in the east coast of northern America, whereas those with decreasing were situated in the west coast of northern America and eastern Asia (Figure 2a). All stomatal models overestimated mean Δ13C values by 1.8‰ to up to 5.0‰ compared to the global TR compilation (Figure 2b), which was higher than the bias reported with the leaf isotopic data set. The amplitude of the trends and the IAV of Δ13C predicted by JULES were lower than those of the observational network (SDs of slopes of 0.02 vs. 0.08‰ year−1, and average IAV of 0.2 ± 0.1‰ vs. 0.5 ± 0.2‰, respectively; Figure 2b,c). It is possible that the uncertainty or variability in f post increased the bias between the TR‐derived ∆13C and the JULES predictions.

FIGURE 2

Tree‐ring‐derived Δ13C trends across the globe over the 1979–2012 period (a). Theil‐Sen's slopes (b) and interannual variations (c) against mean Δ13C values derived from the tree‐ring isotopic measurements and the (summer average) JULES predictions for each of the stomatal models, that is, Jacobs (blue), Leuning (light blue), Medlyn (yellow), and Prentice (red). In (b), the empty and full dots correspond to nonsignificant (p > .05) and significant (p < .05) slopes, respectively

Using the published δ13C data set in bulk leaves and corresponding α‐cellulose TR at the AmeriFlux stations to estimate the post‐photosynthetic fractionation, we found large variations in f post not only across species and sites, but also between the 2 years of available data, with values averaging 4.1 ± 1.1‰ (Table 2). This estimate falls within the higher range of the standard f post values reported in the literature (i.e., 2.1 ± 1.2‰). Using the bulk leaf and stem isotopic data set from the CarboEuropeFlux stations corrected by 1.1‰ to account for difference between bulk and α‐cellulose stem (Guerrieri et al., 2017), we found an apparent f post of around 2.2 ± 1.0‰ (Table 2). While this last estimate is consistent with the standard f post values we used in the ‘global network’ analyses, accounting for the f post calculated at the AmeriFlux stations would result in lower TR‐derived ∆13C than previously estimated, reducing the difference between observed and predicted ∆13C values.

TABLE 2

Apparent post‐photosynthetic fractionations between bulk leaf and α‐cellulose in TR δ13C (f post, ‰) and associated standard deviation (SD, ‰) estimated at the eddy‐covariance sites from the AmeriFlux and CarboEuropeFlux networks

Network	Site	Species	Years	f post	SD
AmeriFlux	US‐Ha1	QURU	2003, 2013	1.97	0.26
		TSCA		6.08	0.25
	US‐Bar	FAGR	2003, 2013	4.32	0.02
		TSCA		5.71	0.68
	US‐DK2	CATO	2002	3.20
		LITU		3.90
	US‐Ho1	PIRU	2003, 2013	5.49	0.15
		TSCA		4.08	1.88
	US‐SP1	PIEL	2013	4.85
		PIPA	2002, 2013	3.64	0.47
	US‐MMS	ACSA	2005	4.36
		LITU		2.74
	US‐Slt	PIEC	2013	4.00
		QUPR		3.00
	Us‐Fuf	PIPO	2014	3.47
	Average			4.05	1.14
CarboEuropeFlux	BE‐Bra	PISY_QURO	2001, 2002	1.88	0.14
	IT‐Col	FAGR	2001, 2002	3.00	0.72
	DE‐Hai	FAGR	2001, 2002	1.85	0.38
	FI‐Hyy	ABAL	2001	2.11
	IT‐Lav	ABAL	2002	3.44
	FR‐LBr	PIPI	2001, 2002	0.72	0.53
	NL‐Loo	PISY	2001, 2002	2.81	0.42
	DK‐Sor	FAGR	2002	0.65
	DE‐Tha	PIAB	2001, 2002	3.43	0.28
	Average			2.20	1.04

We thus tested JULES at the AmeriFlux stations using the TR isotopic records corrected with the individual average f post values reported in Table 2 for each species and site. Consistent with the results of the global network, predictions using the Prentice model better captured the site‐to‐site variations observed in the isotopic record, especially the lowest ∆13C values at the driest site, while also producing average ∆13C values generally closer to the observations (Figures S2 and S3). As the rate of carbon assimilated during photosynthesis depends on leaf c i (Equation 1 and Equations S1–S3) and photosynthesis and transpiration rates are strongly coupled via leaf stomata (Equations 4 and 14), we, therefore, expected stronger predictive skills of the Prentice stomatal scheme to simulate carbon and water fluxes as measured at the top of the canopy. Consistently, we found that the mean values and interannual variations of GPP and ET measured during the growing season at the AmeriFlux and CarboEuropeFlux stations were better reproduced by JULES using the Prentice model (Pearson's r = .5, RMSE = 1.1; Figure S4). Therefore, the Prentice model shows the lowest biases for GPP and ET as well as for the isotopic measurements.

Global variations in Δ13C over 1979–2016 and contributions to these changes

We ran JULES globally using the Prentice stomatal model with the different discrimination models. We found a global average GPP‐weighted Δ13C value for forest ecosystems of 23.8 ± 1.9‰ (mean ± SD) using the simple form of the model, 22.9 ± 1.2‰ incorporating the photorespiratory effect, and 20.9 ± 1.5‰ when additionally incorporating the mesophyll effect. There was, however, large spatial variability in the difference of Δ13C values with or without the photorespiratory or mesophyll effects. The photorespiratory effect increased average Δ13C values in mountainous areas but reduced them in other regions (Figure 3a). The mesophyll effect decreased average Δ13C values everywhere across the globe but the effect was smaller in the tropical forest regions and larger in dry and mountainous areas (Figure 3b). The spatial pattern of mesophyll effect resembled that of ET and T r (Figure S5), with lower impact on Δ13C when (evapo)transpiration was high than when it was low.

FIGURE 3

Differences of average Δ13C values (a) without and with photorespiratory effect (ΔΔ13Cphoto = Δ13Csimple−Δ13Cphoto) and (b) without and with mesophyll effect (ΔΔ13Cmeso = Δ13Cphoto−Δ13Cphoto+meso) over 1979–2016. Average values (c) and trends (d) in Δ13C including both photorespiratory and mesophyll effects over 1979–2016. On the left sides of each panel are the corresponding latitudinal averaged values. In panel (d) only Δ13C trends significant at 90% are shown (p < .10)

The final GPP‐weighted Δ13C values (including both photorespiratory and mesophyll effects) varied strongly across the globe, with higher Δ13C in regions in tropical forests and lower Δ13C in arid and mountainous areas (Figure 3c), consistent with the spatial patterns observed in the leaf isotopic data set (Figure 1c). The trends in Δ13C over the 1979–2016 period also differed across regions. Δ13C increased in the cold high latitudes of the Northern Hemisphere and/or relatively wet areas such as India but decreased in hot and/or dry regions such as central west United States, Patagonia and eastern Russia (Figure 3d), relatively consistent with the trends inferred from the TR isotopic data set (Figure 2a).

Globally, annually Δ13C decreased slightly over the 1979–2016 period when considering the simple discrimination model (−0.0021 ± 0.0004‰ year−1, p < .001; Figure 4a) because of the slight decrease in c i/c a in the global JULES predictions (−0.0001 ± 0.0000 year−1, p < .001; Figure 4b). However, it stayed relatively constant globally when incorporating the photorespiratory effect and additional fractionation due to the diffusion of CO2 within the mesophyll (Figure 4a). The global intrinsic water‐use efficiency (iWUE), which is directly related to c i/c a and thus Δ13C, increased by 0.28 ± 0.01 µmol mol−1 year−1 over the same period (Figure 4c).

FIGURE 4

Global average change in (a) Δ13C for simple version of the discrimination model (dark green), with photorespiratory effect (light blue) and with both photorespiratory and mesophyll effects (dark blue), (b) c i/c a and, (c) iWUE over 1979–2016. The trends and associated standard deviations and P‐value of the trends are reported

When combining the Δ13C trends within different groups of climatic regions (Figure 5a,b), we did not find clear pattern of variation over the globe, except in temperate regions (annual T air of 5–20°C) where Δ13C tended to decrease (up to −0.4‰ over the whole 1979–2016 period), and in cold and wet regions where it increased. In contrast, GPP stayed constant or increased almost everywhere in the globe (Figure S6a), but the rate of increase was slightly lower in cold regions with annual T air lower than 5°C (Figure S6b,c).

FIGURE 5

(a, b) Δ13C trend scores (‰ over the whole 1979–2016 period) and (d, e) Δ13C‐GPP correlation scores for groups of sites with different ranges of annually average T air and β soil values (a, d) or D values (b, e) over 1979–2016. (c) Correlation map between Δ13C and GPP showing only correlations significant at 90% (p < .10). The scores are calculated as the average of the Δ13C trends (a, b) or Δ13C‐GPP correlations (d, e) within each group. The black numbers in the middle of each square correspond to the percentage of data within the group. Only groups with more than 20 grid‐points are considered

Sign, magnitude, and drivers of Δ13C‐GPP relationships over 1979–2016

Large variations in the sign and magnitude of the Δ13C‐GPP relationship were observed between regions (Figure 5c), with negative correlations in tropical Africa or in high northern latitudes such as in Alaska, northern Eurasia and eastern Russia, but positive correlations elsewhere. When combining the Δ13C‐GPP grid‐point correlations within the different groups of climatic regions, we found that Δ13C and GPP were positively related in most environments except in wet‐humid and cold areas (annually averaged T air < 5°C and β soil > 0.6 or D < 0. 8 kPa) and in dry‐warm region (annually averaged T air = 10–15°C and β soil < 0.3) where they tended to be negatively correlated (Figure 5d,e).

The multiple linear regression models showed that Δ13C and GPP values tended to increase with T air but decrease with D in most groups of Δ13C‐GPP correlation, with atmospheric CO2 only slightly contributing to modulate Δ13C and GPP (Table 3). However, while Δ13C was predominantly influenced by D variations, GPP was mainly driven by T air. Around 0.8 ± 0.9% and 55.4 ± 6.3% of the variance in GPP was explained by D and T air, respectively, in group of sites where Δ13C and GPP were negatively related, but only 9.5 ± 0.4% and 25.6 ± 19.5% of GPP variability was explained by the two drivers, respectively, in the other groups. These results suggest that T air had a stronger influence on GPP in cold‐wet and dry‐warm areas. In contrast, 35.2 ± 6.0% and 5.5 ± 0.5% of the variance in Δ13C was explained by D and T air, respectively, in sites with negative Δ13C‐GPP correlations, while 31.8 ± 8.2% and 26.2 ± 4.5% was explained by the two drivers, respectively, in sites with positive correlations.

TABLE 3

Summary statistics for the environmental dependencies of Δ13C (‰) and gross primary production (GPP) (gC m−2 s−1) within groups of Δ13C‐GPP correlations. Standardized fitted coefficients of the linear regression models are reported with percentage of the variance explained by each fixed effect in parenthesis. The coefficient of determination (R 2) is also shown

Correlation group	Climatic range	Variable	CO2	T air	D	R 2
Negative	T air < 5°C and β soil > 0.6	Δ13C	0.01 (<0.01%)	0.25 (5.8%)	−0.63 (35.2%)	.41
		GPP	0.16 (2.7%)	0.89 (62.1%)	0.06 (0.3%)	.65
	T air = 15–20°C and β soil < 0.3	Δ13C	0.01 (<0.01%)	<0.01 (4.8%)	−0.62 (29.2%)	.34
		GPP	0.04 (2.3%)	0.66 (54.4%)	−0.15 (0.3%)	.57
Positive	T air ≥ 5°C or β soil ≤ 0.6	Δ13C	0.01 (<0.01%)	1.42 (23.0%)	−1.51 (26.0%)	.49
		GPP	0.23 (0.4%)	3.05 (39.4%)	−1.47 (9.2%)	.49
Negative	T air < 5°C and D ≤ 0.8 kPa	Δ13C	0.01 (<0.01%)	0.28 (5.8%)	−0.75 (41.2%)	.47
		GPP	0.16 (1.38%)	1.17 (49.7%)	−0.23 (1.9%)	.53
Positive	T air ≥ 5°C or D > 0.8 kPa	Δ13C	0.01 (<0.01%)	1.49 (29.4%)	−1.69 (37.6%)	.67
		GPP	0.26 (0.5%)	1.55 (11.8%)	−1.41 (9.7%)	.22

DISCUSSION

Examining decadal variations in Δ13C of C3 woody plants across the globe and their link to GPP is important for understanding how forest ecosystems have adjusted their physiology with environmental changes over the recent historical period. Previous studies have suggested that (1) Δ13C increased globally at least since 1978 (Keeling et al., 2017) and (2) Δ13C and GPP are positively related to each other (Belmecheri et al., 2014; Peters et al., 2018). Our results based on a new carbon isotopic capability in the land‐surface model JULES challenge these assumptions showing that not only Δ13C trends varies across regions and stayed nearly constant globally over the past decades but that the directionality of the relationship between Δ13C and GPP strongly depends on the climatic regions considered.

Δ13C values and trends vary across regions as a response to environmental conditions

The average Δ13C values predicted by JULES vary depending on the discrimination model considered, with almost 3‰ of difference in Δ13C between the simple and the full models. The large variability in the photorespiratory effect on Δ13C values observed across the globe mainly reflects the impact of atmospheric pressure (and thus elevation) on the photorespiratory term (i.e., −fΓ*/c a, Equation 5), while the one of the mesophyll effect is related to the impact of evaporative demand on the mesophyll term (i.e., −a m(c i − c c)/c a, Equation 5). When exposed to drier conditions, plants reduce both their stomatal (g sc) and mesophyll (g m) conductance (Dewar et al., 2018; Knauer et al., 2020) to minimize water loss during transpiration but at the expense of carbon gain. The reduction of g m in addition to g sc with water stress may reduce Δ13C values in drier environments, resulting in a larger mesophyll effect on Δ13C in dry areas. In contrast, in tropical forests around 20°N and 20°S where the (evapo)transpiration is important, the mesophyll effect is smaller.

Overall, trends in Δ13C varied widely across the globe, with tendencies towards decreasing in temperate and warm regions but increasing in some cold and/or wet areas. Varying temporal patterns of Δ13C trends were also reported in the literature (Andreu‐Hayles et al., 2011; Keller et al., 2017; Lavergne et al., 2017; Martinez‐Sancho et al., 2018), but there were never clearly associated to site environmental conditions. Contrasting Δ13C trends between tree species growing in the same sites have also been observed (Guerrieri et al., 2016; Levesque et al., 2017), which may reflect differences in ecophysiological characteristics between species. These potential differences in Δ13C trends between PFTs were not predicted by JULES because no fixed parameter defined the behavior of forest PFTs in the Prentice stomatal model. Our results show that the sensitivity of Δ13C to rising CO2 is modulated by climate variations and that no single scenario of physiological response to CO2 (e.g., following theoretical framework proposed by Saurer et al. (2004)) can explain this pattern.

Global predicted trend in Δ13C and its contributions differ from Keeling et al. (2017) study

Our simulations using the Prentice model and including discrimination effects from photorespiration and mesophyll conductance suggest that Δ13C stayed relatively constant globally during the 1979–2016 period. Our estimate contrasts with Keeling et al. (2017) findings (referred to as K2017 hereafter), which reported a global increase of Δ13C of 0.014 ± 0.007‰ ppm−1 over 1978–2014 based on atmospheric measurements and a box model. In K2017, the authors assumed a constant Γ*, though it has been shown that Γ* depends on atmospheric pressure and leaf temperature (Bernacchi et al., 2002, 2003). For instance, the increase in global T air of 0.029 ± 0.007°C year−1 over 1978–2014 may have led to an increase of Γ* of 0.055 ± 0.014 ppm year−1 for a range of leaf temperature values of 17–25°C. Such an increase in Γ* counteracted the positive effect of rising CO2 (1.72 ± 0.05 ppm year−1) on Δ13C leading to a lower contribution from the photorespiratory effect to the global Δ13C trend estimated by JULES than by K2017 study (i.e., 0.0015 ± 0.0050‰ ppm−1 vs. 0.0041 ± 0.0014‰ ppm−1).

The authors in K2017 also based their calculations on a linear increase of A of 45% for a doubling of CO2 as suggested by Franks et al. (2013) (equivalent to 6.3% increase for the range of CO2 values from 330 to 377 ppm over 1978–2014), but they assumed constant global mean g m value over time, while g m is expected to decrease with rising atmospheric CO2 (Flexas et al., 2012; Knauer et al., 2020). As a result, the authors found a contribution from the mesophyll effect to the global increasing Δ13C trends of 0.006 ± 0.003‰ per ppm increase of CO2 (using a mesophyll term defined as −(b − a m)/g m A/c a, with A = 9 µmol m−2 s−1 and g m = 0.2 mol m−2 s−1). This estimate strongly differs from that predicted by JULES where the mesophyll effect reduces the global trend in Δ13C by −0.0006 ± 0.001‰ ppm−1. Using a lower g m value at 377 ppm of CO2 than that used by K2017 gives a mesophyll contribution of around 0.002‰ ppm−1 with g m = 0.19 mol m−2 s−1 but of −0.003‰ ppm−1 with g m = 0.18 mol m−2 s−1, more in line with our JULES estimate. Thus, the assumption made in K2017 about the g m value strongly influenced their mesophyll contribution estimate and may have biased their interpretations.

Our modelled and TR‐based estimates of Δ13C trends suggest that changes in Δ13C of C3 plants cannot account for an increase in apparent biospheric discrimination as suggested by K2017. If there were changes to biosphere‐atmosphere isotopic fluxes over the last several decades with an effect on the global δ13CO2 trend that resembled an apparent increase in Δ13C, they must rather derive from non‐photosynthetic processes, such as changes in post‐photosynthetic fractionation, changes in species abundance and distribution or in the proportion of C3:C4 plant productivity, or changes in the residence time of carbon in plants or soils.

The negative Δ13C‐GPP correlation in cold‐wet and tropical African forests reflects the strong influence of T air on photosynthesis in these environments

Overall, Δ13C and GPP tended to be negatively correlated in wet‐cold regions and tropical African forests (where Δ13C remained constant or decreased, while GPP increased over 1979–2016; see also Figures 3d and S6a), but were positively related elsewhere. Because Δ13C and A vary in similar directions with changes in c i/c a and Γ*/c a (see also Text S2), increasing with rising c a and T air but decreasing with rising D, we would expect Δ13C and canopy‐scale A (i.e., GPP) to be positively correlated everywhere. Given that T air and D are positively related to each other (Pearson's r = .70, p < .001), it is possible that the dominant influences of T air on GPP (via A) and of D on Δ13C in particular in cold and wet‐humid regions both contributed to the decoupling between Δ13C and GPP. The stronger impact of T air on GPP than on Δ13C could be due to the additional effect of rising T air on the maximum carboxylation capacity (V cmax) or electron transport (J max), which do not impact Δ13C. It is also possible that it reflects the positive influence of increasing T air on the leaf area (Luyssaert et al., 2007; Piao et al., 2006) or on the growing season length which promotes carbon uptake in those regions (Cai & Prentice, 2020; Piao et al., 2009).

The spatial distribution of the Δ13C‐GPP correlations found in European forests (i.e., negative north of 60°N but positive elsewhere) is consistent with Shestakova et al. (2019) results in the region using Δ13C and TR growth. The authors of this study also found a change from negative to positive correlations in the Scandinavian area over the course of the 20th century very likely due to increasing drying. Their results suggest that if the atmospheric drying observed over the past century (Grossiord et al., 2020) was to continue, it is possible that the Δ13C‐GPP relationship in relatively wet and cold regions may switch from negative to positive because of the increasing influence of D on GPP. Thus, studying variations in Δ13C and GPP together may help documenting the areas where the drivers of photosynthesis could change in the coming decades due to climate change.

Limitations of the study

We acknowledge several limitations in our approach to examining Δ13C mean values and trends and Δ13C‐GPP relationships. First, the Δ13C values across the globe predicted by JULES (including both photorespiratory and mesophyll effects) were around 0.8‰ higher in average than those derived from the global leaf network, and up to 1.6‰ higher than those from the TR network. In addition, the magnitudes of the Δ13C trends were almost four time lower in the predictions than in the TR isotopic network and the IAV in Δ13C was largely underestimated by the model. The apparent model bias in Δ13C values may result from uncertainties in both the predicted and observation‐derived Δ13C.

On the one hand, it is possible that the ‘least‐cost’ Prentice model overestimate c i (or c c) because it ignores the effect of soil drying on stomata aperture, leading to an overestimation of Δ13C. A recent study showed that the ratio of c i to c a could be reduced by up to 2% globally when considering the impact of soil water stress on the stomatal activities (Lavergne, Sandoval, et al., 2020). These findings suggest that c i (and c c) should decrease with a reduction of plant‐available soil water, resulting in lower predicted Δ13C values in particular in dry areas, more in line with the observed values. The use of yearly average (instead of monthly) c a values to predict Δ13C may also have led to an overestimation of the simulated Δ13C values because of the strong seasonal cycle of c a (i.e., lower c a values in summer than in winter). Given that the seasonal variations in c a have increased since the 1950s in particular in the northern high latitudes (Graven et al., 2013) but were not accounted for here, we could expect the IAV in Δ13C to be lower in the predictions than in the TR observations. Finally, we assumed a constant g m/g s following recent studies (Flexas & Carriqui, 2020; Yiotis & McElwain, 2019) to estimate c c (Equation (1), (9a), (9b)) while it is still unclear whether g m responds differently or more strongly than g s to environmental changes. For instance, (all else equal) if g m/g s increases (e.g., g m changes at a higher rate than g s), c c (and thus c c/c a) should decrease, resulting in an even more negative effect of g m on Δ13C trend. In contrast if g m/g s decreases (e.g., g s changes at a higher rate than g m), c c (and thus c c/c a) should increase, reducing the negative effect of g m on Δ13C trend. Thus, our assumption of a constant g m/g s may have affected the predicted Δ13C values. Nevertheless, a recent study found no significant difference in g m/g s values across long‐term drought‐stress treatments (Ma et al., 2021), increasing confidence about our assumption.

On the other hand, biases in the observation‐derived Δ13C may have increased the discrepancies between predicted and observed Δ13C values. Since the yearly average δ13CO2 values used to calculate Δ13C from the leaf and TR isotopic observations are lower than those during the growing season, the data‐derived Δ13C values may have been overestimated at least in the Northern Hemisphere where the δ13CO2 seasonal cycle varies between 0.2‰ in tropical areas to up to 0.8‰ in high latitudes (Scripps CO2 program). Also, we assumed a constant f post to correct the TR isotopic timeseries, while our analyses based on isotopic data from the AmeriFlux stations suggest that f post varies strongly across sites, species and even years, and may be in average twice higher than the standard f post values often used in the literature (4.1‰ vs. 2.1‰). Using the higher f post estimates would result in higher TR‐derived ∆13C values, reducing the difference between observed and predicted ∆13C values, and thus the apparent bias in the predictions.

Despite all the above‐mentioned limitations, the spatial patterns of Δ13C mean values predicted by JULES were in good agreement with the leaf isotopic observations and with Cornwell et al. (2018) study based on a data‐driven statistical approach, suggesting that the model captures relatively well the spatial variability in Δ13C values across the globe. Since the average Δ13C values inferred from the JULES simulations for the TR sites decreased over 1979–2016 (−0.006‰ year−1, p < .001), whereas those estimated from the observed TR Δ13C timeseries stayed nearly constant (p = .097), it is, however, possible that JULES has a negative bias in the global trend in Δ13C. One reason for this bias may be that the slight decrease in the predicted c i/c a values over 1979–2016 (reduction of 0.5%, p < .001) is not realistic and that c i/c a should rather stay constant over time as suggested by other studies (Frank et al., 2015; Wong et al., 1979). Nevertheless, JULES predicted an almost proportional increase of iWUE compared to atmospheric CO2 over 1979–2016 (21% vs. 20%, respectively), consistent with the rate of increase estimated by recent studies for similar periods (Adams et al., 2020; Lavergne et al., 2019; Mathias & Thomas, 2021; Soh et al., 2019), suggesting that despite biases to predict IAV and trend in Δ13C, the carbon uptake per unit of water lost is relatively well constrained in JULES.

CONCLUSIONS

This is the first study investigating decadal changes in Δ13C and documenting the links between Δ13C and GPP across the globe using a data‐model approach. Δ13C predicted by JULES increases in some cold and wet environments but decreases in temperate regions over the 1979–2016 period. These patterns of variations show the strong coupling between the water and carbon cycles across the globe. Globally, predicted Δ13C stayed nearly constant over the studied period, which differs from the large secular increase in Δ13C reported in K2017 based on atmospheric measurements. The difference between K2017 and this study mainly lies on the contrasting simulated contributions from photorespiratory and mesophyll effects to the global trend. Our results suggest that studying Δ13C and GPP together may inform about the most influential drivers of photosynthesis and, thus, help understanding temporal changes in the land carbon uptake.

AUTHOR CONTRIBUTIONS

A.L. designed the research, analyzed the data, and wrote the first draft. D.H. and R.G. provided the stable carbon isotope data from leaves and TRs in the CarboEuropeFlux and AmeriFlux stations, respectively. R.J.O. helped in the setup and run of JULES on the NERC JASMIN platform using the Rose/Cylc suite control system. I.C.P. and H.G. provided additional suggestions on the research plans. All co‐authors contributed to improve the manuscript.

CONFLICT OF INTEREST

All authors declare that they have no conflicts of interest.

CODE AVAILABILITY STATEMENT

The model code and the files needed for running it are available from the Met Office Science Repository Service (MOSRS; https://code.metoffice.gov.uk/trac/jules/; registration required). The results presented in this paper were obtained by running JULES vn5.6 branch with a new carbon isotopic modelling capability (code available after registration with MOSRS at https://code.metoffice.gov.uk/trac/jules/browser/main/branches/dev/alienorlavergne/vn5.6_jules_Cisotopes). The runs were performed with the Rose suite u‐bx886 (https://code.metoffice.gov.uk/trac/roses‐u/browser/b/x/8/8/6). The Python code developed to create the main manuscript figures is available at https://github.com/Alielav/GCB_Lavergneetal_2021.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.