Low sensitivity of gross primary production to elevated CO₂ in a mature eucalypt woodland

The response of mature forest ecosystems to a rising atmospheric carbon dioxide concentration (C_a) is a major uncertainty in projecting the future trajectory of the Earth's climate. Although leaf-level net photosynthesis is typically stimulated by exposure to elevated C_a (eC_a), it is unclear how this stimulation translates into carbon cycle responses at the ecosystem scale. Here we estimate a key component of the carbon cycle, the gross primary productivity (GPP), of a mature native eucalypt forest exposed to free-air CO₂ enrichment (the EucFACE experiment). In this experiment, light-saturated leaf photosynthesis increased by 19 % in response to a 38 % increase in C_a. We used the process-based forest canopy model, MAESPA, to upscale these leaf-level measurements of photosynthesis with canopy structure to estimate the GPP and its response to eC_a. We assessed the direct impact of eC_a, as well as the indirect effect of photosynthetic acclimation to eC_a and variability among treatment plots using different model scenarios.

At the canopy scale, MAESPA estimated a GPP of 1574 g C m⁻² yr⁻¹ under ambient conditions across 4 years and a direct increase in the GPP of +11 % in response to eC_a. The smaller canopy-scale response simulated by the model, as compared with the leaf-level response, could be attributed to the prevalence of RuBP regeneration limitation of leaf photosynthesis within the canopy. Photosynthetic acclimation reduced this estimated response to 10 %. After taking the baseline variability in the leaf area index across plots in account, we estimated a field GPP response to eC_a of 6 % with a 95 % confidence interval (−2 %, 14 %). These findings highlight that the GPP response of mature forests to eC_a is likely to be considerably lower than the response of light-saturated leaf photosynthesis. Our results provide an important context for interpreting the eC_a responses of other components of the ecosystem carbon cycle.

1 Introduction

Forests represent the largest long-term terrestrial carbon storage (Bonan, 2008; Pan et al., 2011). The atmospheric carbon dioxide concentration (C_a) has increased significantly since the beginning of the industrial era (Joos and Spahni, 2008), but the increase would have been considerably larger without forest carbon sequestration, which is estimated to have offset 25 %–33 % of recent anthropogenic CO₂ emissions (Le Quéré et al., 2018). C_a is projected to continue to increase by 1–5 µmol mol⁻¹ yr⁻¹ into the future (IPCC, 2014), but the rate of this rise depends on the magnitude of the forest feedback on C_a. At the leaf scale, the direct physiological effects of rising C_a are well understood: elevated C_a (eC_a) stimulates plant photosynthesis (Kimball et al., 1993; Ellsworth et al., 2012) and reduces stomatal conductance (Morison, 1985; Saxe et al., 1998), which together increase leaf water-use efficiency (De Kauwe et al., 2014). These physiological responses at the leaf scale can increase ecosystem carbon uptake, which in turn may result in increased carbon storage in the ecosystem, mitigating the rise in C_a. However, projecting the response of the terrestrial carbon sink to future increases in C_a is a major uncertainty in models (Friedlingstein et al., 2014), highlighting an urgent need to make greater use of data from manipulative experiments at the leaf scale to inform terrestrial biosphere models (Medlyn et al., 2015).

Our understanding of ecosystem responses to eC_a relies on both experiments and observations. However, results from different types of studies show some important areas of disagreement (Fatichi et al., 2019). At the global scale, satellite data provide evidence of a strong greening trend over the last 20 years, indicating an increase in leaf area and/or above-ground biomass, which has been attributed to the gradual increase in CO₂ (Donohue et al., 2009, 2013; Yang et al., 2016; Zhu et al., 2016). A positive response of carbon uptake/greenness is also found in manipulative eC_a open-top chamber experiments with young trees (Eamus and Jarvis, 1989; Curtis and Wang 1998; Saxe et al., 1998; Medlyn et al., 1999) and ecosystem-scale FACE experiments in young, aggrading forest stands (Ainsworth and Long, 2005; Norby et al., 2005; Ellsworth et al., 2012; Walker et al., 2019). In contrast, individual-tree experiments with mature trees (> 30 years old) have found relatively small responses of tree growth to eC_a despite an apparent increase in leaf photosynthesis (Dawes et al., 2011; Sigurdsson et al., 2013; Klein et al., 2016). Moreover, tree-ring studies indicate an apparent lack of stimulation of vegetation growth in mature forests over the last century (Peñuelas et al., 2011; Silva and Anand, 2013; van der Sleen et al., 2014). These studies raise important questions about how mature ecosystems will respond to eC_a.

The Eucalyptus FACE experiment (EucFACE; Australia) is the first replicated, ecosystem-scale experiment where a mature native forest has been experimentally subjected to eC_a and provides a valuable case study to assess the response of a mature forest to eC_a under field conditions (Ellsworth et al., 2017). Results from the first 5 years (2013–2018) of leaf gas exchange measurements showed a consistent stimulation of leaf-level light-saturated net photosynthesis (A) of 19 % (Ellsworth et al., 2017; Wujeska-Klause et al., 2019). Nevertheless, the increase in A did not lead to a detectable change in above-ground growth (Ellsworth et al., 2017). These experimental results are consistent with empirical evidence arising from tree-ring studies (Peñuelas et al., 2011; Silva and Anand, 2013; van der Sleen et al., 2014) and also with experimental evidence from individual mature trees (Körner et al., 2005; Dawes et al., 2011; Klein et al., 2016).

As a first step towards reconciling the eC_a responses of leaf photosynthesis and above-ground growth in this experiment, we quantify how the whole canopy carbon uptake, or gross primary productivity (GPP), increased under eC_a. The response of the GPP is important because it provides a reference point against which to compare the response of other components of the ecosystem carbon balance, such as above-ground growth. It needs to be quantified explicitly because the response of the GPP to eC_a may be quite different to that of the leaf net photosynthesis. The leaf-level response of photosynthesis to eC_a is usually measured on sunlit leaves under saturating light (Ainsworth and Rogers, 2007). As a result, these leaf-level eC_a responses largely reflect the responses of the photosynthesis rate when limited by maximum Rubisco activity (V_cmax). However, depending on the canopy architecture and the ambient light conditions, the canopy could have many shaded leaves, which would mean that the emergent rate of photosynthesis could actually be limited by RuBP regeneration (J). RuBP-regeneration-limited photosynthesis has a smaller response to eC_a than Rubisco-limited photosynthesis (Ainsworth and Rogers, 2007), resulting in a smaller response of the GPP than the leaf photosynthesis under saturating light.

The transition from RuBP-regeneration-limited to Rubisco-limited photosynthesis of the canopy is determined by the ratio of the maximum capacities for RuBP regeneration and Rubisco activity, J_max and V_cmax (Friend, 2001; Zaehle et al., 2014; Rogers et al., 2017). Wullschleger (1993) reported a J_max:V_cmax ratio of 2, which has been widely adopted in models (e.g. Wang et al., 1998; Luo et al., 2001; Rogers et al., 2017). However, recent studies have suggested that the J_max:V_cmax ratio varies systematically across forest ecosystems and can range from 1 to 3 (Kattge and Knorr, 2007; Ellsworth et al., 2012; Kumarathunge et al., 2018). A lower J_max:V_cmax ratio results in more frequent RuBP regeneration limitation of photosynthesis, which reduces the response of the GPP to eC_a.

It is difficult to directly measure the eC_a effect on the GPP. In some previous eC_a experiments, the GPP has been estimated by scaling up from leaf-level measurements using a canopy model. Wang et al. (1998) and Luo et al. (2001) both used the tree array model, MAESPA, which can simulate the radiative transfer within and between tree crowns and can be parameterized to describe the spatial locations and sizes of trees in eC_a experiments. In these previous applications of MAESPA, the direct response of GPP to eC_a was consistently half of that observed at the leaf level due to the large contribution of RuBP-regeneration-limited photosynthesis to the GPP (Wang et al., 1998; Luo et al., 2001). However, the direct effect of eC_a on photosynthesis was modified by two major indirect effects: (1) when LAI increased under eC_a, the additional leaf area amplified the GPP response by up to 60 %; and (2) the downregulation of photosynthesis under eC_a, or photosynthetic acclimation (Long et al., 2004; Ainsworth and Rogers, 2007; Rogers, et al., 2017). Under long-term exposure to eC_a, some plants have been observed to reduce nitrogen allocation to Rubisco, which results in a decrease in the photosynthetic capacity (Gunderson and Wullschleger, 1994). The average decrease of V_cmax among plants in FACE experiments was found to be 13 % for all species and 6 % for trees (Ainsworth and Long, 2005). Both Wang et al. (1998) and Luo et al. (2001) tested the impact of photosynthetic acclimation and showed a moderate reduction in the canopy GPP (5 %–6 %) due to photosynthetic acclimation (10 %–20 %).

Following Wang et al. (1998) and Luo et al. (2001), we used MAESPA (Duursma and Medlyn, 2012) to estimate the canopy GPP at EucFACE in ambient and elevated C_a treatments. The model has previously been evaluated using leaf- and tree-scale measurements from EucFACE (Yang et al., 2019). Here, we first parameterized the model using physiological and structural data measured during the experiment. Then, we quantified the response of the canopy GPP to eC_a and partitioned this response into the direct stimulation of the GPP and the indirect effects of the photosynthetic acclimation and variation of LAI. The overall goal of this study was to estimate the magnitude of the response of the forest canopy GPP to eC_a in order to provide a baseline against which to compare changes in other components of the ecosystem carbon balance.

2 Methods

2.1 Site

The EucFACE experiment (technical details in Gimeno et al., 2015) is located in western Sydney, Australia (33.617^∘ S, 150.741^∘ E). It consists of six circular plots, each of which has a diameter of 25 m, enclosing 15–25 mature forest trees (referred to as “rings” hereafter). The rings are divided into two groups: control (with ambient C_a; 390–400 µmol mol⁻¹ during the study period) and experimental (eC_a; +150 µmol mol⁻¹). The tree canopy is dominated by Eucalyptus tereticornis Sm. which are ∼20 m in height and have a basal area of ∼24 m² ha⁻¹. The site receives a mean annual precipitation of 800 mm yr⁻¹, a mean annual photosynthetically active radiation (PAR) of 2600 MJ m⁻² yr⁻¹, and has a mean annual temperature of 17 ^∘C.

2.2 Model

The MAESPA model is a process-based tree-array model (Wang and Jarvis, 1990) that calculates canopy carbon and water exchange (https://bitbucket.org/remkoduursma/maespa/src/Yang_et_al_2019/, last access: 4 July 2019). At each 30 min time step, the model simulates the radiative transfer, photosynthesis, and transpiration of individual trees mechanistically. The soil moisture balance can be calculated dynamically, but here we chose to improve accuracy by using soil moisture as an input to the model (Duursma and Medlyn, 2012).

The model represents the tree canopy as an array of tree crowns. The location and dimensions of each crown are specified based on on-site measurements (see Sect. 2.3.2 below). Calculations of carbon and water fluxes are made for each tree crown, which is divided into six layers. Here, it was assumed that crowns are represented by an ellipsoidal shape and that leaf area is uniformly distributed across layers within the tree crown. The leaf angles were assumed to follow a spherical distribution to ensure consistency with the method used to estimate leaf area index (LAI) in Duursma et al. (2016). Within each layer, the model evaluates the radiation transfer and leaf gas exchange at 12 grid points such that each crown is represented by a total of 72 grid points. The radiation intercepted at each grid point is calculated for direct and diffuse components by considering shading from the upper crown and surrounding trees, solar angle (zenith and azimuth), and light source (diffuse or direct). Penetration by direct radiation to each grid point is used to estimate the sunlit and shaded leaf area at each grid point. The radiation intercepted by the fraction of sunlit and shade foliage is then used to calculate the leaf gas exchange.

The gas exchange sub-model combines the leaf photosynthesis model of Farquhar et al. (1980) with a stomatal optimization model, following Medlyn et al. (2011). Stomatal conductance is modelled as follows:

$$\begin{array}{}\text{(1)}& g_{\mathrm{s}}=\mathrm{1.6}\cdot (\mathrm{1}+{\displaystyle \frac{g_{\mathrm{1}}}{\sqrt{D}}})\cdot {\displaystyle \frac{A_{\mathrm{net}}}{C_{\mathrm{a}}}},\end{array}$$

where g_s is the stomatal conductance to water vapour (mol m⁻² s⁻¹), g₁ is a parameter that represents the g_s sensitivity to photosynthesis (kPa^0.5; see definition in Medlyn et al., 2011), A_net is the net CO₂ assimilation rate (µmol m⁻² s⁻¹), C_a is the atmospheric CO₂ concentration (µmol mol⁻¹), and D is the vapour pressure deficit (kPa). The factor 1.6 converts the conductance of CO₂ to that of H₂O.

The impact of soil moisture on g_s is represented via an empirical function that links soil water availability to g₁ following (Drake et al., 2017):

$$\begin{array}{}\text{(2)}& g_{\mathrm{1}}=g_{\mathrm{1},\max }({\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\min }}{{\mathit{\theta }}_{\max }-{\mathit{\theta }}_{\min }}}{)}^q,\end{array}$$

where the g_1,max is the maximum g₁ value; θ is volumetric soil water content (%); θ_max and θ_min are the upper and lower limit within which θ has impact on g₁, respectively; and q describes the nonlinearity of the curve. The equations used to calculate A_net are given in the Supplement (Supplement S1, Eqs. S1–S6).

Following Yang et al. (2019), MAESPA considers a non-stomatal limitation to the biochemical parameters J_max and V_cmax at high D:

$$\begin{array}{}\text{(3)}& V_{\max }=V_{\max ,t}(\mathrm{1}-c_D\cdot D),\end{array}$$

where V_max,t is the J_max or V_cmax at a given leaf temperature (Supplement S1), and c_D is a fitted parameter (Table 1). This relationship is empirical and fitted to data collected in EucFACE. Incorporating this relationship was shown to improve the predicted photosynthesis by the leaf gas exchange model (Yang et al., 2019).

Combining Eqs. (1)–(3) and (S1)–(S6) yields the g_s and A_net of each grid point, which is then multiplied by the leaf area at each grid point and summed to give the whole-tree photosynthesis. Photosynthesis of individual trees is then summed to give the whole-canopy photosynthesis.

2.3 Model parameterization

2.3.1 Meteorological forcing

The model is driven by in situ PAR, wind speed, air temperature, D, and soil moisture measurements from 2013 to 2016 (Figs. 1, 2). Each ring has a set of PAR (LI-190, LI-COR, Lincoln, NE, USA), wind speed (WINDCAP Ultrasonic WMT700 Vaisala, Vantaa, Finland), humidity, and temperature sensors (HUMICAP^® HMP 155 Vaisala, Vantaa, Finland) in the centre of the ring above the canopy at 23.5 m. The PAR, air temperature, and relative humidity were measured every 5 min in each ring and were then gap-filled by linear interpolation and aggregated to 30 min mean time slices across all six rings (Fig. 1). D was calculated from temperature and humidity measurements.

Figure 1 Meteorological data measured at the site during the 2013–2016 period. Panels show (a) the daily mean vapour pressure deficit (D) with the shaded area marking the maximum and minimum of the day, (b) the daily mean air temperature (T_air) with the shaded area marking the maximum and minimum of the day, (c) the daily maximum photosynthetically active radiation (PAR), and (d) the monthly total precipitation. Note that precipitation has no direct impact in the model but modifies stomatal conductance via the change in soil moisture.

Two levels of C_a were used in the model according to the measured C_a (LI-840, LI-COR, Lincoln, NE, USA). The ambient C_a was gap-filled (there were < 10 d in total during the 4 years, with gaps due to a power outage) and aggregated to 30 min mean time slices from the 5 min measurements across the three ambient rings (rings 2, 3, and 6). The eC_a was processed in the same way but using data from the experimental rings (rings 1, 4, and 5).

The volumetric soil water content (θ) was used as an estimate of plant water availability and was taken every 20 d using neutron measurements at 25 cm intervals (503DR Hydroprobe, Instroteck, NC, USA) and averaged to the top 150 cm (Fig. 2). There were two probes in each ring and the average of these probes was used to represent the ring average for each measurement date. θ was updated on the days of measurements and thus not gap-filled.

Figure 2 (a) The leaf area index (LAI) and (b) soil volumetric water content (θ) used to drive the model. The LAI was estimated in each ring from measurements of understorey PAR and smoothed using a generalized additive model following Duursma et al. (2016). θ was measured using neutron probes in the top 150 cm biweekly (Gimeno et al., 2018). Each line colour indicates a different plot. Red shows elevated CO₂ plots (treatment), whereas blue shows ambient CO₂ plots (control). The x axis ticks mark the start of each year.

2.3.2 Canopy structure

Trees in MAESPA were represented by their actual location, height, and crown size to mimic the realistic effects of shading. Tree location, crown height, crown base, and stem diameter were measured in January 2013 at the start of the experiment. For each ring, a time series of the LAI was obtained based on measurements of above- and below-canopy PAR (Duursma et al., 2016). This LAI represents the plant area index, which includes the woody component as well as leaves and does not account for clumping. In order to retrieve the actual LAI, we assumed a constant branch and stem cover (0.8 m² m⁻²) based on the lowest LAI during November 2013 when the canopy shed almost all leaves. Thus, the LAI used in this study was the plant area index estimate from Duursma et al. (2016), which was less 0.8 m² m⁻² (Fig. 2a). As the LAI is the only parameter besides the soil moisture that differed by ring, the canopy structure (i.e. the LAI and its distribution) was the major driver of inter-ring variability.

The total leaf area (m²) of each ring was calculated as the product of the LAI and the ground area of each plot (491 m²). This total leaf area (LA) was then assigned to each tree based on an allometric relationship between the total leaf area (m²) and the diameter at breast height (DBH; cm). The allometric relationship was derived from data in the BAAD database (Falster et al., 2015) for Eucalyptus trees grown in natural conditions with DBH values less than 1 m to match the characteristics of EucFACE. In total, this database yielded a total of 66 observations with which to estimate the relationship between the LA and DBH:

$$\begin{array}{}\text{(4)}& L_{\mathrm{allom}}=a\cdot {\mathrm{DBH}}^b,\end{array}$$

where L_allom is the theoretical leaf area based on an allometric relationship with the DBH. The values obtained via fitting for a and b were 492.6 and 1.8 respectively, with a root mean square error of 14.4 (m²) and a correlation coefficient of 0.83. Equation (4) was used to assign the total LA of each ring to each tree using the following steps: (i) the L_allom for each tree was calculated based on the DBH; (ii) the L_allom was summed to obtain a total LA for each ring; and (iii) the fractional contribution of each tree to the ring total LA was calculated. The total LA based on the LAI was then assigned to each tree based on this fraction.

The crown radius was calculated using a linear function with the DBH based on measurements made in August 2016. The data consisted of the DBH and crown radius (one on a north–south axis and one on an east–west axis) of four trees in each ring. The crown radius measurements were averaged by tree and used to fit a linear model with the DBH. The estimated slope and intercept of the relationship are 0.095 (m cm⁻¹) and 0.765 (m) respectively.

MAESPA also considered the shading from surrounding trees outside the rings. However, no measurements of locations or diameters were available for the trees surrounding the rings. Therefore, a total of 80 surrounding trees were arbitrarily assumed to form two uniform and circular layers around each ring. They were assigned the mean height, mean crown radius, and mean leaf area estimated from all trees in EucFACE. Except for shading, the surrounding trees have no impact on the trees within the rings. Ring 1 is shown in Fig. S1 in the Supplement as an example of the representation of canopy structure in MAESPA.

2.3.3 Physiology

The physiological parameters were estimated from field gas exchange measurements as described below. The data were collected with portable photosynthesis systems (Li-6400, LI-COR, Inc., USA). The only parameter found to differ between ambient and elevated C_a rings was V_cmax,25 (V_cmax at 25 ^∘C; Ellsworth et al., 2017). Hence, all other parameters (e.g. the temperature responses of photosynthesis and respiration) were estimated by combining all data across CO₂ treatments. Fitted parameter values are given in Table 1.

A set of temperature-controlled photosynthesis–CO₂ response (A−C_i) curves was measured at different leaf temperatures (20–40 ^∘C) under saturating light in February 2016. The data set was used to quantify the temperature dependencies of J_max and V_cmax by fitting a peaked Arrhenius function (Eq. S5) to the measurements. We assumed that these temperature response functions applied throughout the period of the study.

Light- and temperature-controlled A−C_i curves were also measured in the morning for 10 field campaigns from 2013 to 2016. All A−C_i curves were started at the growth C_a of 395 µmol mol⁻¹ or 545 µmol mol⁻¹ (depending on the eC_a treatment) with a saturating light of 1800 µmol m⁻² s⁻¹ and a flow rate of 500 µmol s⁻¹ with temperature controlled to a constant based on the seasonal temperature. These data were used to estimate J_max and V_cmax at 25 ^∘C using the “fitaci” function in the “plantecophys” R package (Duursma, 2015), utilizing the measured temperature responses of J_max and V_cmax described in the previous paragraph to correct to 25 ^∘C.

Repeated gas exchange measurements were made on the same leaves in the morning and afternoon under prevailing field conditions and saturating light (photon flux density = 1800 µmol m⁻² s⁻¹) on four occasions in 2013 (“diurnal”; Gimeno et al., 2018). To expand the diurnal data set, we obtained the points from A−C_i curves at field C_a and combined the two data sets. These data were used to estimate the g₁ parameter in the stomatal conductance model (Eq. 1) using the “fitBB” function in the plantecophys R package (Duursma, 2015). One g₁ value was fitted to the data from each treatment and date. The g₁ values were then regressed against θ measured in each treatment group to estimate the impact of soil moisture availability on leaf gas exchange, following Eq. (2). The g₁ values were related to the nearest measurements of θ (within 2 weeks). There was no rainfall between g₁ and θ measurement dates. Equation (2) was fitted to this data set using the non-linear least squares method (Fig. 3).

Figure 3 The impact of the soil moisture content (θ) in the top 150 cm on stomatal regulation. The g₁ parameter values are fitted to measurements of leaf gas exchange grouped by month and treatment. Red dots are fitted to data from elevated rings while blue are ambient rings. Error bars indicate the standard errors of the fitted values. The grey line shows the fit of Eq. (2) to the data.

The dark respiration rate of foliage, R_dark, was measured at least 3 h after sunset at a range of leaf temperatures (14–60 ^∘C) in February 2016 (also with the Li-6400). The temperature dependence of R_dark was fitted to all of the measured data with Eq. (S6) using the non-linear least squared method. Light responses of photosynthesis were measured on two trees from each ring in October 2014 (Crous et al., unpublished). This data set was used to constrain the light response parameters (α_J and θ_J) in Eq. (S4). Details of fitting the light response curves are provided in the Supplement (Supplement S1).

2.4 Model simulations and analysis

MAESPA was used to simulate radiation interception and gas exchange of all six rings between 1 January 2013 and 31 December 2016 on a 30 min basis. The model simulated the 30 min gross primary production (GPP) of each tree, which was then summed for all trees in each ring to get the total annual GPP for each ring and year.

Four different sets of simulations were used to estimate carbon uptake under ambient and eC_a conditions and to identify the key limiting factors on the canopy GPP response to eC_a. Firstly, we carried out a simulation of leaf-scale (“leaf scenario”) photosynthesis with measured meteorological data but fixed physiological data (g₁=3.3 kPa^0.5, $V_{\mathrm{cmax},\mathrm{25}}=\mathrm{91}$, and $J_{\max ,\mathrm{25}}=\mathrm{159}$ µmol m⁻² s⁻¹). This simulation aimed to quantify the CO₂ response of Rubisco-limited and RuBP-limited photosynthesis at the leaf scale. This calculation was made using the “Photosyn” function in the plantecophys R package (Duursma, 2015). This function implements the leaf gas exchange routine used in MAESPA.

Secondly, MAESPA was run for all six rings with ambient C_a and with V_cmax,25 from ambient measurements (“ambient scenario”). The results of this simulation were used to calculate the GPP of each ring under ambient conditions. The ambient GPP values were also used to evaluate the inherent variability among the rings.

Thirdly, all six rings were simulated with eC_a and V_cmax,25 based on measurements from ambient rings (“elevated scenario”). The results of this simulation were compared to those from the ambient scenario to illustrate the instantaneous response of canopy GPP to eC_a in each ring and year. This simulation also quantifies the variation of the GPP response to eC_a across rings and years.

Lastly, we simulated the response of the three rings exposed to eC_a (rings 1, 4, and 5) using the V_cmax,25 and eC_a measured from these elevated rings (“field scenario”). Results from the field scenario were used for two analyses: (i) to compare the GPP from the field scenario to that of the three rings from the elevated scenario (i.e. eC_a and ambient V_cmax,25), which allows us to quantify the impact of photosynthetic acclimation (i.e. due to a reduction in V_cmax); and (ii) to calculate the difference in the GPP between the three ambient rings in ambient scenario and elevated rings in the field scenario to estimate the response of the GPP to eC_a in the field.

Table 1 Summary table of parameter definitions, units, and sources used in this study.

3 Results

Figure 4 summarizes the results from measurements and the different simulations conducted in this study. It demonstrates that the impact of eC_a diminishes as calculations are scaled from the instantaneous leaf-level response (A_inst) to the long-term canopy response (GPP_field) and the various feedback effects are accounted for. Each row in Fig. 4 is explained in detail in the following paragraphs.

Figure 4 The response of photosynthesis to eC_a on different scales and limited by different factors. In summary, from top to bottom, the figure demonstrates how a large increase in leaf photosynthesis can diminish into a statistically nonsignificant change in the canopy GPP under eC_a. Entries from top to bottom are as follows: (a) A_inst is the instantaneous response of leaf photosynthesis to eC_a obtained from A-C_i measurements in ambient rings (error bars indicate 95 % CI); (b) A_c is the modelled response of Rubisco-limited leaf photosynthesis, assuming no downregulation, averaged over the range of diurnal air temperatures experienced during the experimental period; (c) A_J is the modelled response of RuBP-regeneration-limited leaf photosynthesis; (d) GPP_inst is the direct effect of eC_a on the canopy GPP, modelled with MAESPA, assuming no downregulation of photosynthesis and averaged across all six rings; (e) A_long is the long-term response of leaf photosynthesis to eC_a obtained from leaf photosynthesis measured at treatment CO₂ concentrations (see Ellsworth et al., 2017) – this value is different from A_inst because it incorporates photosynthetic acclimation; (f) GPP_long is the effect of eC_a on the canopy GPP once the measured downregulation of V_cmax is taken into account; (g) LAI is the measured difference in the average LAI between eC_a and ambient C_a rings over the experiment period (data from Duursma et al., 2016); and (h) GPP_field is the GPP response modelled with MAESPA comparing the three elevated rings with the three ambient rings. The bars represent model outputs, and points represent observations. See the text for further explanation.

3.1 Instantaneous C_a response of photosynthesis at the leaf and canopy scale

The mean instantaneous C_a response of leaf-level photosynthesis (A_inst) was +33 % (Fig. 4a). This response ratio was calculated from ∼600 light- and temperature-controlled A−C_i curves measured in the ambient rings. From the curves, we extracted the photosynthesis at 400 and 550 C_a (µmol mol ⁻¹) and calculated the instantaneous C_a effect as their ratio. This approach allows for an estimation of the direct CO₂ response independent of the impact of photosynthetic acclimation.

By contrast, the modelled direct GPP response to eC_a was considerably less (just +11 %) as shown in Fig. 4d (“GPP_inst”). This canopy response rate was calculated by comparing the modelled GPP of all six rings under ambient and elevated C_a conditions (“ambient” vs. “elevated” scenario). As a result, this direct canopy GPP response also excludes the impact of photosynthetic acclimation.

Our results show that the major reason for the difference between the direct leaf and canopy photosynthesis responses to eC_a is the relative contributions from Rubisco-limited and RuBP-regeneration-limited photosynthesis (cf. Fig. 4b and c). Figure 5 shows that the response of photosynthesis to eC_a is considerably higher when Rubisco activity limits photosynthesis (A_c) than when RuBP regeneration limits photosynthesis (A_J). When averaged over the range of leaf temperatures experienced during the 4 years of the experiment, the A_c response to eC_a on average (+26 %; Fig. 4b) is larger than that of A_J (+10 %; Fig. 4c). Leaf gas exchange measurements were taken in saturating light (1800 µmol m⁻² s⁻¹) and, thus, are mostly Rubisco-limited. Therefore, the observed response rate of A_inst is close to that of A_c.

Figure 5 The modelled C_a response of Rubisco-limited leaf photosynthesis (A_c) and RuBP-regeneration-limited leaf photosynthesis (A_J) to leaf temperature (T_leaf). The responses are calculated for temperatures during the 2013–2016 period. Parameters are as given in Table 1, except that V_cmax,25 and g₁ were assumed to be constant for clarity ($g_{\mathrm{1}}=\mathrm{3}.$ kPa^0.5 and $V_{\mathrm{cmax},\mathrm{25}}=\mathrm{90}$ µmol m⁻² s⁻¹).

At the canopy scale, a large fraction of the modelled canopy photosynthesis is limited by RuBP regeneration. In Fig. 6, we show the distribution of A_c and A_J during the 4 years of simulation as calculated by MAESPA. On average, 70 % of the canopy photosynthesis is limited by RuBP regeneration under ambient conditions (“ambient scenario”). The high fraction of A_J is partly a consequence of the relatively low ratio of J_max,25 to V_cmax,25 (J:V ratio) which was estimated to be 1.7 (Table 1). In Fig. 7, we estimated the PAR level at which Rubisco activity becomes limiting to leaf photosynthesis. The transition point from Rubisco-limited to RuBP-regeneration-limited photosynthesis was calculated from the leaf gas exchange sub-model by assuming a constant C_a (390 µmol mol⁻¹), D (1.5 kPa), g₁ (3.3 kPa^0.5), and V_cmax,25 (90 µmol m⁻² s⁻¹), but varying leaf temperature. As shown, under these conditions, when the temperature is 25 ^∘C and the J:V ratio is 1.7, Rubisco activity limits photosynthesis only when incident PAR is greater than 1800 µmol m⁻² s⁻¹. Using a higher J:V ratio such as the commonly used value of 2 would decrease the saturating PAR value at which photosynthesis becomes Rubisco-limited. We ran additional simulations assuming a J:V ratio of 2 and found that, with this ratio, MAESPA estimated 48 % of photosynthesis to be RuBP-regeneration-limited under ambient conditions and a direct GPP response of 15 % (data not shown).

Figure 6 Distribution of average annual photosynthesis limited by Rubisco activity and RuBP regeneration in bins of PAR (30 µmol m⁻² s⁻¹). As calculated by MAESPA across all rings from 2013 to 2016. The histogram was constructed by calculating the photosynthesis (either limited by Rubisco or RuBP) falling into each bin for every 30 min in the “ambient scenario”. These values were then summed to each year and ring and averaged over six rings and 4 years.

Figure 7 Estimated PAR value at which limitation to photosynthesis shifts from RuBP regeneration to Rubisco at different leaf temperatures and J:V ratios. Rubisco limitation occurs at PAR values above the curves; RuBP regeneration limitation occurs below the curves. The curves were calculated using the Photosyn function in the plantecophys R package (Duursma, 2015). Parameters other than PAR and T_leaf were assumed to be constant: C_a of 390 µmol mol⁻¹, D of 1.5 kPa, g₁ of 3.3 kPa^0.5, and V_cmax,25 of 90 µmol m⁻² s⁻¹. The temperature and light dependences of photosynthesis were assumed to be the same as in MAESPA. The grey line was predicted by assuming $J_{\max ,\mathrm{25}}=\mathrm{153}$ µmol m⁻² s⁻¹ (i.e. a J:V ratio of 1.7). This J:V ratio was observed consistently in EucFACE across campaigns and rings. The red line was predicted by assuming $J_{\max ,\mathrm{25}}=\mathrm{180}$ µmol m⁻² s⁻¹ (i.e. a J:V ratio of 2). This J:V ratio has been commonly reported and used in other studies. The horizontal dashed line shows the PAR = 1800 µmol m⁻² s⁻¹ at which leaf-level measurements of EucFACE were made. Note the log scale of the y axis. The dashed curves are based on quantum yield of electron transport (α_J; mol mol⁻¹) and the convexity of light response of RuBP, θ_J, unitless values from the CABLE model (Haverd et al., 2018).

The shape of the light response curve also determines the transition point from RuBP-limited to Rubisco-limited photosynthesis. We explored this effect by investigating the effect of varying the convexity, θ_J, which is assumed to be the same as the convexity of overall photosynthesis. At EucFACE, we estimated this parameter to be 0.48 from light-response curves of photosynthesis collected on site, indicating a shallow curvature and a high light saturation point, in contrast to the more commonly assumed 0.85, representing a steeper curvature and a lower light saturation point. Using a value of 0.85 for θ_J resulted in a much lower PAR being required for photosynthesis to became Rubisco-limited (dashed curves in Fig. 7). With a θ_J of 0.85 and a J:V ratio of 1.7, MAESPA estimated 40 % of photosynthesis to be RuBP-regeneration-limited under ambient conditions and a direct GPP response of 16 % (data not shown). With a θ_J of 0.85 and a J:V ratio of 2, MAESPA estimated just 34 % of photosynthesis to be RuBP-regeneration-limited under ambient conditions and a direct GPP response of 18 % (Fig. S2). Thus, the simulated CO₂ response of canopy carbon uptake depends heavily on the parameterization of light response and the J:V ratio.

3.2 Acclimation of photosynthesis

The above calculations are made considering only the instantaneous response of photosynthesis to eC_a. However, photosynthetic acclimation has been observed at the leaf scale (Ellsworth et al., 2019), and will also reduce the response of the GPP to eC_a at the canopy scale. At the leaf-level, photosynthesis measured in the elevated rings after 5 years of treatment (A_long) was 19 % higher than that measured in ambient rings (Fig. 4e; Ellsworth et al., 2017). Thus, A_long accounts for the photosynthetic acclimation in the elevated rings after 4 years of exposure to eC_a. A_long is considerably smaller than A_inst (19 % vs. 33 %; Fig. 4a and e), indicating a large effect of photosynthetic acclimation on the eC_a response of light-saturated photosynthesis.

Accounting for the impact of photosynthetic acclimation in MAESPA, by using the V_cmax from elevated rings (“field” vs. “ambient” scenarios), reduced the response of the GPP to C_a from 11 % to 10 % (GPP_long; Fig. 4f). As such, the photosynthetic acclimation had a relatively modest impact on the modelled annual GPP in the model. The small impact of photosynthetic acclimation on canopy photosynthesis relative to the effect on leaf photosynthesis can be explained by the fact that the leaf photosynthesis data are measured under saturating light and, thus, are typically Rubisco-limited, so a reduction in V_cmax had a large effect. In contrast, at the canopy scale, much of the photosynthesis was limited by RuBP regeneration and was largely unaffected by a reduction in V_cmax.

3.3 Influence of LAI

The realized GPP response to eC_a also depends on the canopy structure, specifically the LAI. In this experiment, there was no significant change in the LAI with eC_a ($-\mathrm{4}\mathit{\%}\pm \mathrm{5}\mathit{\%}$; Fig. 4g; see also Duursma et al., 2016). The effect of eC_a on the LAI was calculated as the average effect between the elevated and ambient annual mean LAI. However, there was inherent variability in the LAI across the rings (Fig. 2a), which does not fundamentally change the effect of eC_a but requires a detailed analysis of the potential effects of natural variability on the response to eC_a.

The small pretreatment difference in the LAI across rings gives rise to a range of estimates for the GPP response to eC_a in the field (6 %±8 %; Fig. 4h). This result is explored further in Fig. 8, which combines the results from “ambient”, “elevated”, and “field” scenarios. The average GPP across all six rings under ambient C_a was 1574 g C m⁻² yr⁻¹ over the 4-year simulation (“ambient scenario”; Fig. 8). However, there was significant variability in the ambient GPP across rings, related in part to the inherent variability in the LAI across rings. We characterized the pre-existing differences in the LAI by the initial LAI (LAI_i), measured on 26 October 2012. These initial values are low, as they are measured immediately before the seasonal leaf flush, but characterize the difference in the LAI across rings over the full experimental period. Rings 1 and 4 (both experimental rings) have the lowest LAI_i (< 0.3 m² m⁻²) and thus the lowest average GPP under ambient conditions (1206 g C m⁻² yr⁻¹). Ring 5 (the other experimental ring) has the second highest LAI_i (∼0.4 m² m⁻²) and also the highest GPP under ambient conditions (2359 g C m⁻² yr⁻¹). Therefore, the variability among rings in the ambient GPP (SD of 15 %) is larger than the modelled direct effect of C_a on GPP, which is similar in all rings (+11 %).

Figure 8 The 4-year average GPP of all six rings under ambient and eC_a conditions plotted against the initial leaf area index (LAI_i). The LAI_i is the LAI measurement taken on the 26 October 2012 and is a proxy for the inherent variation among the rings. For all six rings, the estimated GPP is shown for ambient C_a (blue) and eC_a (red). Crosses indicate the GPP from simulations in which the C_a is varied, and squares indicate the GPP under field conditions. The flat bars on the right-hand side of the plot indicate the average ambient C_a GPP for ambient rings only (the average of blue squares) and the average eC_a GPP for elevated rings only (the average of red squares). Dashed lines indicate the average ambient C_a (the average of blue crosses) and the eC_a GPP across all six rings (the average of red crosses). Thus, the flat bars mark the modelled response without inter-ring variability, whereas the dashed lines mark the modelled realized response, including inter-ring variability.

Owing to the variability among rings represented by the LAI_i, the estimated mean GPP response to eC_a across the experimental rings has a sizable confidence interval (±8 %, Fig. 4h). The actual eC_a response was estimated as an average effect between the ambient and elevated GPP values considering the impacts of photosynthetic acclimation and inter-ring variability. The average GPP of experimental rings under field conditions (eC_a) was estimated to be 1698 g C m⁻² yr⁻¹, whereas the average GPP of control rings under field conditions (ambient C_a) was 1599 g C m⁻² yr⁻¹, which is an increase of 6 % as shown in the Fig. 4h. The variation of the annual average GPP of the control and experimental groups (blue and red squares in Fig. 8) are consequently represented by the CI in Fig. 4h.

4 Discussion

We have shown how a large response of leaf-level photosynthesis to eC_a diminishes when integrated to the canopy scale, according to the synthesis of 4 years of leaf measurements at EucFACE with a stand-scale model, MAESPA. We estimated that the canopy GPP of a mature Eucalyptus woodland under ambient C_a conditions varied from 1084 to 2129 g C m⁻² yr⁻¹ by ring and year with a mean of 1574 g C m⁻² yr⁻¹. The model, constrained by site measurements, predicted that once scaled to the canopy, the response of the GPP to eC_a only increased by 6 % (95 % CI of ±8 %) compared with the 19 % (95 % CI of ±5 %) observed in leaf-scale measurements. We were able to quantify the response of the GPP to eC_a and attribute the reduction in the response to various factors including (i) Rubisco limitations vs. RuBP regeneration limitations to photosynthesis, (ii) photosynthetic acclimation, and (iii) inter-ring variability in the LAI. Together, these findings provide valuable insights into the relative importance of each factor and help close a key knowledge gap in our understanding of how mature forests respond to eC_a.

4.1 Performance of MAESPA under ambient conditions

The ambient GPP of EucFACE estimated by MAESPA was comparable to that measured with eddy covariance in similar evergreen eucalypt forests in southeastern Australia. At a nearby eddy covariance site (< 1 km), Renchon et al. (2018) estimated the ecosystem GPP from eddy covariance to be 1561 g C m⁻² yr⁻¹ from 2013 to 2016; this is within the range estimated for the ambient rings in this study, although this latter site and the EucFACE site are not the same in terms of canopy structure and LAI. Furthermore, our version of MAESPA was evaluated against leaf photosynthesis and whole-tree sap flow measurements in EucFACE (R² of 0.77 and 0.8 respectively; Yang et al., 2019). These comparisons indicate that MAESPA is a useful tool to explore the canopy carbon uptake, and the predicted GPP could provide a baseline for future studies.

4.2 RuBP-regeneration-limited photosynthesis

Our results show that the canopy GPP at EucFACE was predominantly limited by RuBP regeneration. The reason for the frequent RuBP regeneration limitation is that the measured J:V ratio was relatively small in EucFACE (1.7), and stomata tend to close at midday when light levels are higher and Rubisco limitation is expected (Gimeno et al., 2015). A lower J:V ratio increases the PAR threshold required for the photosynthesis model to switch between the RuBP regeneration limitation and the Rubisco limitation (from < 1000 to < 1800 µmol m⁻² s⁻¹; Fig. 7). Previous studies have highlighted the need to consider the J:V ratio for a correct prediction of the CO₂ response (Long et al., 2004; Zaehle et al., 2014; Rogers et al., 2017). However, as shown by Zaehle et al. (2014), Medlyn et al. (2015), and Rogers et al. (2017), current models differ with respect to their predictions of the transition from RuBP-regeneration-limited to Rubisco-limited photosynthesis, suggesting the uncertainty of the predicted CO₂ response of the GPP could be reduced by using a realistic J:V ratio.

Previous modelling studies applying MAESPA to eC_a experiments both assumed a higher J:V ratio (2) and estimated a higher GPP response to eC_a, presumably due to less frequent RuBP regeneration limitation (Wang et al., 1998; Luo et al., 2001). A J:V ratio of 2 was suggested by Wullschleger (1993) and has been used in many modelling studies (e.g. the seven terrestrial biosphere models assessed by Rogers et al., 2017 all assumed a J:V ratio of between 1.9 and 2). Global terrestrial biosphere models such as JULES and others frequently estimate J_max on the basis of this ratio (e.g. Clark et al., 2011). However, the relatively low J:V ratio observed at EucFACE is not unique. At the Duke Forest FACE site in the US, Ellsworth et al. (2012) reported a J:V ratio of ∼1.7, which is the same as that estimated for EucFACE. Kattge and Knorr (2007) analysed V_cmax and J_max values from 36 species across the world and found a low J:V ratio (< 1.8) in herbaceous, coniferous, and broad-leaved species. Most recently, Kumarathunge et al. (2018) studied the variation in the J:V ratio in data sets obtained from around the globe and found that it declined with increasing growing season temperature. The ratio varied from 2.5 in tundra environments to less than 1.5 in tropical environments. The value of 1.7 observed at EucFACE falls within this prediction for the prevailing growth temperature at this site. Thus, the inclusion of this relationship between the J:V ratio and temperature will be important for capturing the GPP response to eC_a globally.

We also found that the convexity of the light response of photosynthesis affected the predicted GPP response to eC_a (Fig. 7). The parameter value we fitted to data measured in situ (θ_J=0.48) is lower than the value commonly assumed in the models (e.g. 0.7 in Bonan et al., 2011). Note that some model studies assume that θ_J is lower than the convexity of overall photosynthesis (typically over 0.8; e.g. 0.9 in Medlyn et al., 2002; 0.85 in Haverd et al., 2018). Here, we assumed that the convexity of the electron transport rate and overall photosynthesis are the same (see Supplement S1 for details). Nonetheless, our relatively low θ_J value (< 0.7) is not unique, as it is also supported by a number of studies on different species around the world (Ögren, 1993; Valladares et al., 1997; Lewis et al., 2000; Hjelm and Ögren, 2004). The inclusion of a higher θ_J value would predict a much higher direct GPP response to eC_a (e.g. 16 % vs. 11 % in this study), as higher θ_J results in a large proportion of the GPP being Rubisco-limited. This finding calls for careful examination of the light response of photosynthesis, which has a large effect on the predicted eC_a response

4.3 Photosynthetic acclimation

Some degree of photosynthetic acclimation (i.e. a long-term reduction of V_cmax under eC_a) has been widely reported in FACE studies and has been attributed to a reduction in the leaf nitrogen concentration (Saxe et al., 1998; Ainsworth and Long, 2005). The response of the GPP to eC_a would be linearly related to V_cmax if photosynthesis were mostly limited by Rubisco activity. Photosynthetic acclimation was responsible for the reduced response of leaf-scale light-saturated photosynthesis from 33 % (A_inst) to 19 % (A_long). However, this reduction in V_cmax only translated into a ∼2 % reduction in the GPP modelled by MAESPA. Wang et al. (1998) also showed that photosynthetic acclimation (−21 % in V_cmax) reduced the modelled canopy GPP by only 6 % due to RuBP regeneration being the primary limitation of canopy photosynthesis. Thus, these findings suggest that photosynthetic acclimation may only have a small effect on the GPP response to eC_a when canopy photosynthesis is mostly RuBP-regeneration-limited. Therefore, this response is consistent with the hypothesis that the reduction in V_cmax represents a reallocation of nitrogen to optimize the nitrogen use efficiency under eC_a (Chen et al., 1993; Medlyn et al., 1996).

4.4 Constraining the carbon balance response to eC_a

At EucFACE, after 4 years of eC_a treatment, there was no evidence of increased above-ground tree growth (Ellsworth et al., 2017). Furthermore, the trees at EucFACE also did not show any significant change in the LAI (Duursma et al., 2016). The relatively small response of the GPP and the effect of ring-to-ring variation provides important context for these statistically nonsignificant responses of tree growth at the stand scale at EucFACE. Firstly, the effect size calculated for the GPP of +11 % (+169 g C m⁻² yr⁻¹) constrains the likely effect size for plant growth and other components of the ecosystem carbon balance and is a more useful baseline for comparison than the response of light-saturated leaf photosynthesis (+19 % = 299 g C).

Secondly, the inherent ring-to-ring variation in this natural forest stand is larger than the GPP response, which highlights the importance of considering both the effect size and variability in the observations rather than just focusing on statistical significance. It is important to note that the EucFACE site could be considered relatively homogeneous for a mature woodland. The site is flat, trees appear to be similar with respect to age, and almost the entire overstorey is comprised of a single species. In addition, plots were carefully sited to minimize variation in basal area. However, there are small-scale variations in soil type, depth, and nutrient availability that cause variation in the LAI. This scale of variation is likely to present in other natural forests, and, indeed, other studies on mature trees also note that background variability can contribute to a lack of statistically significant findings (Fatichi and Leuzinger, 2013; Sigurdsson et al., 2013). We highlight the need to focus on effect size and its uncertainty, rather than the dichotomous significant/nonsignificant approach when evaluating experimental results from native forests.

4.5 Implications for terrestrial biosphere models

Seven terrestrial biosphere models (TBMs) were used to predict the GPP and LAI responses to eC_a in advance of the EucFACE experiment (Medlyn et al., 2016). The predicted eC_a responses of the GPP ranged from +2 % to +24 % across the seven models, while the predicted responses of the LAI ranged from +1 % to +20 %. With our results, it is possible to disprove some of the assumptions made in these model simulations and identify directions for model improvement. The model with the lowest GPP response (CLM4-P) assumed very strong downregulation of photosynthesis owing to phosphorus limitation. However, this downregulation was not observed here. The models with the highest GPP responses (GDAY, O-CN, and SDGVM) had a J:V ratio of 2, which is higher than that observed at EucFACE, and also had a positive feedback to the GPP via increased LAI (+5 % to 15 %), which did not occur (Duursma et al., 2016). The model rendering the most similar prediction for the GPP response to eC_a to the output of MAESPA incorporating empirical observations was the CABLE model. This latter model predicted an eC_a response of GPP of ∼12 % with a large proportion of RuBP-regeneration-limited photosynthesis, both of which are similar to the findings in this study. Future TBMs may benefit from incorporating a more realistic representation of the relative contribution of RuBP-regeneration-limited and Rubisco-limited photosynthesis to the GPP. For instance, adding the temperature dependency of the J:V ratio could help capture the variation in the J:V ratio globally (e.g. Kumarathunge et al., 2018).

Our study provides a number of process-based insights that can be used to improve model performance both qualitatively and quantitatively. Our modelling exercise is also a major contribution to the understanding of the EucFACE experiment as it quantifies the amount of extra carbon input into the system by canopy-level photosynthesis and, thus, provides a reference for assessing the impacts of eC_a on growth and soil respiration. Finally, our study highlights that the eC_a effect on canopy-scale GPP may be considerably lower than the effect on the photosynthesis of the light-saturated leaves, due to the contrasting relative limitations to photosynthesis operation and the different scales. In future work, our GPP estimates will be used as an input to calculate the overall effect of eC_a on the carbon balance at the whole EucFACE site.