Assessing terrestrial biogeochemical feedbacks in a strategically geoengineered climate

Geoengineering by injecting sulfur dioxide (SO2) into the lower stratosphere has been suggested to reduce anthropogenically induced warming. While impacts of such geoengineering on climate have been investigated in recent decades, few modeling studies have considered biogeochemical feedbacks resulting from such intervention. This study comprehensively characterizes responses and feedbacks of terrestrial ecosystems, from an ensemble of coupled high-resolution Earth system model climate change simulations, under the highest standard greenhouse gas scenario with an extreme geoengineering mitigation strategy. Under this strategy, temperature increases beyond 2020 levels due to elevated anthropogenic carbon dioxide (CO2) were completely offset by the SO2 injection. Carbon cycle feedbacks can alter the trajectory of atmospheric CO2 levels by storing or releasing additional carbon on land and in the ocean, thus moderating or amplifying climate change. We assess terrestrial biogeochemical feedbacks to climate in response to geoengineering, using model output from the Stratospheric Aerosol Geoengineering Large Ensemble (GLENS) project. Results indicate terrestrial ecosystems become a stronger carbon sink globally because of lower ecosystem respiration and diminished disturbance effects under geoengineering. An additional 79 Pg C would be stored on land by the end of the twenty-first century, yielding as much as a 4% reduction in atmospheric CO2 mole fraction without marine biogeochemical feedbacks, compared to the high greenhouse gas scenario without geoengineering.


Introduction
Rising global mean surface temperature along with increasing anthropogenic greenhouse gas emissions have been observed since the last century (IPCC 2014), and future projections of increasing global mean surface temperature remain even if anthropogenic emissions are reduced (Steffen et al 2018). To prevent continued warming that could cause devastating damage to natural ecosystems and human socio-economic activities (Hoegh-Guldberg et al 2018), various climate intervention strategies, such as solar radiation management (SRM), have been proposed to offset the risks of warming (Crutzen 2006, Shepherd 2009). SRM includes various techniquessurface albedo approaches (e.g. brightening of human settlements), cloud-albedo enhancement (e.g. marine cloud brightening), stratospheric aerosols (e.g. sulfate aerosols oxidized from gaseous hydrogen sulfide or SO 2 ), and space-based techniques (e.g. space reflectors placed in low Earth orbit)-that reduce incoming solar radiation by intentionally increasing Earth's surface albedo (Shepherd 2009). Among the SRM techniques, stratospheric aerosol proposals draw attention widely because of their relatively low cost for mitigating the surface temperature warming and their evenly-distributed cooling capability around the world (Solar Radiation Management Governance Initiative 2011). Inspired by the surface cooling effects of large volcanic eruptions that expel gases into the stratosphere (Robock 2000), most stratospheric aerosol research in past decades focused on stratospheric SO 2 injections (e.g. Robock et al 2008, Kravitz et al 2013, Tjiputra et al 2016, Xia et al 2016 since the gas can be released by large volcanic eruptions into the stratosphere and react to form sulphate aerosols, which can be circulated through the stratospheric winds and dim the incoming insolation. Extensive analyses of Earth System Model (ESM) simulations have addressed climate impacts of stratospheric SO 2 injection, including stratospheric ozone depletion (Tilmes et al 2008), weakening of monsoons (Robock et al 2008, Tilmes et al 2013, acid deposition and ocean acidification (Kravitz et al 2009, Moreno-Cruz andKeith 2013), repartitioning of direct and diffuse radiation (Gu et al 2003, Mercado et al 2009, and suppressed precipitation (Bala et al 2008). Terrestrial biogeochemical (BGC) responses have received limited attention in such simulations, particularly for optimized aerosol injection strategies or for high greenhouse gas forcing scenarios using fully coupled ESMs (Keller et al 2014, Tjiputra et al 2016, Xia et al 2016, Cao and Jiang 2017, Keith et al 2017, Dagon and Schrag 2019, Duan et al 2020. A recent geoengineering modeling activity known as the Stratospheric Aerosol Geoengineering Large Ensemble (GLENS) project  simulated the climate response to stratospheric SO 2 injections designed to maintain three temperature goals at their 2020 levels, under the Representative Concentration Pathway 8.5 (RCP8.5) emission scenario (Riahi et al 2011): (1) the global mean surface temperature, (2) the interhemispheric surface temperature gradient, and (3) the equator-to-pole surface temperature gradient. The simulations were conducted with the Community Earth System Model (CESM) version 1 with the Whole Atmosphere Community Climate Model (WACCM) as its atmospheric component (Mills et al 2017), the Community Land Model version 4.5 with prescribed time-evolving distributions of vegetation consistent with RCP8.5 and interactive carbon and nitrogen cycles as the land component (Oleson et al 2013), the Los Alamos Sea Ice Model (Community Ice CodE version 4) as the sea ice component (Hunke and Lipscomb 2008), and the Parallel Ocean Program version 2 as the ocean component (Danabasoglu et al 2012). SO 2 was injected in the lower stratosphere at optimized locations to achieve the three temperature goals (Kravitz et al 2017) from 2020 to 2099 and to avoid overcooling of the tropics and undercooling of the poles that lead to continued Arctic summer sea ice loss (Moore et al 2014, Tilmes et al 2016). The large number of ensemble members, a longer simulation period, higher model grid resolutions, effects of sulfate aerosols on stratospheric chemistry, and a high anthropogenic emission scenario through the coupled CESM simulations enable GLENS to serve the needs of a more comprehensive examination of geoengineering impacts, including effects on atmospheric chemistry and land biogeochemical feedbacks to the Earth system. While the fully coupled GLENS simulation results have been used to investigate a wide variety of physical climate responses to aerosol geoengineering (e.g. Tilmes et al 2018, Fasullo et al 2018, Cheng et al 2019, the responses and feedbacks of terrestrial ecosystems to this geoengineered climate have not been evaluated. This study complements those as the first to comprehensively characterize the responses and feedbacks of terrestrial ecosystems to rapidly rising CO 2 levels, producing strong CO 2 fertilization responses, without significant changes in global mean surface temperature or consequent warming-induced changes in precipitation after 2020. Our results not only quantify the terrestrial BGC feedbacks in GLENS but also illustrate the importance of terrestrial ecosystems on future climate change in a strategically geoengineered climate.  (Kravitz et al 2017) starting at year 2020. The SO 2 injection rate was adjusted annually during 2020-2099, using a feedback algorithm that maintained global mean surface temperature, the interhemispheric surface temperature gradient, and the equator-to-pole surface temperature gradient at their 2020 levels (Macmartin et al 2014, Kravitz et al 2016Kravitz et al , 2017. Because only the first three (of 20) ensemble members in the baseline experiment were extended to at least year 2097, the first three ensemble members from the baseline experiment during 2010-2019 (called 'BASE' hereafter) and during 2020-2097 (called 'CTRL' hereafter), and the first three ensemble members from the geoengineering experiment during 2020-2097 (called 'GEOENG' hereafter) are analyzed to fairly compare the differences in climate and terrestrial biogeochemical feedbacks between the baseline and geoengineering experiments, as well as those differences with respect to the present climate. In this study, we apply the two-tailed Student's t-test to evaluate if changes between two experiments are significant (p < 0.1). In performing the t-test, we properly adjusted for lag-one autocorrelation by computing an equivalent sample size that effectively decreases the number of degrees of freedom in determining the standard error of the mean for each of the time series (see supplementary information section 2 for details (available online at stacks.iop.org/ERL/15/104043/mmedia)). Reported correlation coefficients for two variables from spatial maps of differences between experiments are evaluated using the linear Pearson's correlation without removing the global mean (called uncentered), following the method described in IPCC (2001).

Ecoregions
We define 13 ecoregions in this study by clustering climate, soil, and topography information into regions that were fitted to the International Geosphere-Biosphere Programme (IGBP) ecoregion definitions (Townshend 1992

Adjusted trajectories of atmospheric CO 2 and surface temperatures
The carbon budget in the terrestrial ecosystem can be expressed, according to the Intergovernmental Panel on Climate Change (Watson et al 2000) as where GPP is the gross primary production, NPP the net primary production, R a the autotrophic respiration, NEP the net ecosystem production, R h the heterotrophic respiration, NBP the net biome production, and disturbance includes anthropogenic emissions due to land cover and land use changes such as fires and crop harvest. In this study, GPP, NPP, NEP, NBP, R a , and R h are direct outputs from GLENS while disturbance is estimated by subtracting NBP from NEP. Positive values for GPP, NPP, NEP, and NBP indicate carbon gains in land while negative values denote carbon losses to the atmosphere; contrarily, positive values for R a and R h represent carbon losses to the atmosphere and larger negative values imply more carbon remaining in land. We utilize NBP to examine the carbon storage in terrestrial ecosystems. A positive NBP value indicates atmospheric CO 2 is stored in the terrestrial ecosystem while a negative NBP value represents carbon releases from the terrestrial ecosystem to the atmosphere. A fixed ratio of 1:2.13 (Clark 1982, O'Hara 1990) is used in this study to convert the amount of carbon (unit: Pg C) released from or sequestered in the terrestrial ecosystem to the equivalent amount of airborne CO 2 mole fractions (unit: ppm). Simulations in GLENS were driven by the atmospheric CO 2 mole fraction specified in the RCP8.5 emission scenario (Meinshausen et al 2011, Riahi et al 2011). The simulations were conducted with an active global land carbon cycle; however, the terrestrial BGC feedbacks were not incorporated into the coupled modeling system. In addition, the marine BGC feedbacks were excluded in the GLENS experiments. We assume the NBP variations are consistent with the atmospheric CO 2 for the CTRL experiment. Hence, differences in simulated NBP between GEOENG and CTRL are attributed to ecosystem responses due to SO 2 injections in the lower stratosphere. A positive value of these differences, i.e. a larger accumulated NBP in GEOENG than in CTRL indicates more carbon is stored from the atmosphere to land, thus the atmospheric CO 2 level is lowered; on the contrary, the atmospheric CO 2 level is higher when less carbon is sequestered in land due to lower accumulated NBP in GEOENG than in CTRL. To construct an adjusted atmospheric CO 2 trajectory that accounts for this feedback, we compute the adjusted CO 2 airborne mole fraction every year according to the equation (2) where t is the time step, C the adjusted atmospheric CO 2 airborne mole fractions, F the atmospheric CO 2 airborne mole fractions obtained from GLENS outputs, and B the terrestrial BGC feedbacks from experiment GEOENG and CTRL. The changes in the atmospheric CO 2 trajectory alter the radiative forcing, resulting in a different surface temperature trajectory. The changes in surface temperatures due to atmospheric CO 2 adjustments are approximated through an impulse response function tuned to the mean of CMIP5 models (Boucher andReddy 2008, Hoffman et al 2014). We use this function to compute CO 2induced radiative forcing changes and thus the associated temperature trajectory adjustments. Differences between the simulated global mean surface temperature in GLENS and the adjusted one are compared to estimate a reduced amount of SO 2 injection, with an injection rate of 10 Tg SO 2 per year for 1 K cooling .

Global climate changes
Global mean surface temperatures are projected higher in CTRL than in BASE due to rising greenhouse gas levels as previously reported (IPCC 2014) (figure 2(a)). In GEOENG, slightly higher mean surface temperatures in the Sahara Desert and at midto-high latitudes compared to BASE (figure 2(b)) are compensated by lower surface temperatures in the central US, parts of South America, Central Asia, southern India, and eastern Australia. The global mean surface temperatures in GEOENG are lower compared to CTRL despite having the same greenhouse gas trajectories as a result of stratospheric aerosol injections (figure 2(c)). Higher precipitation rates are projected in most regions in CTRL than in BASE due to increased water vapor content in the lower troposphere and atmospheric circulation, both of which are induced by the increased temperature in the CTRL simulation (Collins et al 2013); lower precipitation rates are simulated in the Amazon Basin, Chile and Argentina, Uruguay, western Congo, southern Africa, some southern parts of Europe, and Indonesia  (j)). In GEOENG, cooler temperatures, reduced evapotranspiration, and aerosol-cloud interactions (suppressed precipitation due to the aerosol indirect effect (Albrecht 1989)) lead to a drier climate compared to BASE over the southern Sahel and South Africa, India, Southeast Asia, and parts of the boreal zone across Eurasia and northeastern North America (figure 2(e)). Global mean precipitation rates in GEOENG are generally smaller than that in CTRL (figure 2(f)) due to lower temperatures and induced suppression of precipitation, with the exception of the southern Amazon In terms of downward solar radiation at the surface, direct radiation is 109.4 ± 0.3 W m −2 for BASE in 2019, 111.5 ± 1.0 W m −2 for CTRL and 79.8 ± 1.2 W m −2 for GEOENG in 2097; diffuse radiation is 61.6 ± 0.1 W m −2 for BASE in 2019, 58.3 ± 0.4 W m −2 for CTRL and 83.0 ± 0.8 W m −2 for GEOENG in 2097. More direct radiation is projected in CTRL compared to BASE in most regions except for India, northern and central Africa, central South America, and the southeast US (figure 2(g)).
Reductions in cloudiness (figure S1) are the primary explanation for such direct radiation changes (spatial correlation = −0.76), which is consistent with the results from Yu et al (2016). Nevertheless, regions like northeast Australia experience increased downward direct radiation despite increasing cloudiness due to enhanced precipitation (figure 2(d)) that removes suspended aerosols from the atmosphere (i.e. reduced Figure 2. Changes of spatial distributions between CTRL and BASE (top row), between GEOENG and BASE (middle row), and between GEOENG and CTRL (bottom row) for surface temperature (K, (a)-(c)), precipitation (mm day −1 , (d)-(f)), total downward direct solar radiation at the surface (W m −2 , (g)-(i)), and total downward diffuse solar radiation at the surface (W m −2 , (j)-(l)). The spatial distribution of CTRL is from the 2020-2097 time-averaged results without geoengineering while that of GEOENG is from the 2020-2097 time-averaged results with geoengineering applied. The spatial distribution of BASE is from the 2010-2019 time-average results. Light grey stippling (dots) indicates regions where the change is significant using the Student's t-test (p < 0.1). aerosol optical depth (AOD), figure S2). Reduced diffuse radiation at the surface (figure 2(j)) is associated with the reduction of total cloudiness in most regions. On the other hand, GEOENG shows decreased direct radiation and increased diffuse radiation at the surface with respect to BASE (figures 2(h) and (k)) because of increased amounts of aerosols suspended in the atmosphere due to SO 2 injections. An unexpected increase in surface radiation in the dust region of central Australia in GEOENG compared to BASE (figure 2(h)), similar to that described above when comparing CTRL with BASE, is caused by lower dust emission (figures S3 and S4) because of wetter climate conditions (figure 2(e)). In this region, the spatial correlation is −0.94 between precipitation and coarse mode dust AOD, and that between precipitation and Aitken mode dust AOD is −0.59. Differences in downward solar radiation at the surface between GEOENG and CTRL are associated with higher aerosol amounts due to SO 2 injections, which block direct radiation and enhance diffuse radiation (figures 2(i) and (l)).

Global terrestrial biogeochemical responses
Photosynthesis rates (FPSN) and GPP for both CTRL and GEOENG substantially increase in most of the world because of rising atmospheric CO 2 levels with respect to BASE (figures 3(a), (b), (d), and (e)). At high latitudes, FPSN and GPP are lower in BASE and in GEOENG than in CTRL because higher surface temperatures in CTRL (figure 2(c)) thaw permafrost regions, lengthen growing seasons and enable enhanced photosynthesis (figures 3(a) and (c)). At mid and low latitudes, lower precipitation and less diffuse radiation in the Amazon Basin and southern Africa produce lower FPSN and GPP in CTRL than in BASE (figures 3(a) and (d)). In addition, lower FPSN and GPP in India and central Africa, where land cover is dominated by croplands and savannas, are subjected to lower precipitation in BASE and in GEOENG than in CTRL (figures 2(d) and (f) , different spatial patterns are found between FPSN and GPP across the Amazon Basin for changes between GEOENG and CTRL (figures 3(c) and (f)). Since the terrestrial nitrogen cycle has been suggested to be a critical factor for the carbon cycle response (Tjiputra et al 2016), we investigate the spatial pattern of net nitrogen mineralization (supplementary figure S6) and find that the inconsistent spatial patterns between FPSN and GPP are where net nitrogen mineralization is smaller in GEOENG than in CTRL. This is mainly because Figure 3. Changes of spatial distributions between CTRL and BASE (top row), between GEOENG and BASE (middle row), and between GEOENG and CTRL (bottom row) for photosynthesis rates (µmol m −2 s −1 , (a)-(c)), gross primary production (kg C m −2 yr −1 , (d)-(f)), net primary production (kg C m −2 yr −1 , (g)-(i)), and net biome production (kg C m −2 yr −1 , (j)-(l)). The spatial distribution of CTRL is from the 2020-2097 time-averaged results without geoengineering while that of GEOENG is from the 2020-2097 time-averaged results with geoengineering applied. The spatial distribution of BASE is from the 2010-2019 time-average results. Light grey stippling (dots) indicates regions where the change is significant using the Student's t-test (p < 0.1). lower temperatures in GEOENG produce lower net nitrogen mineralization rates, reducing conversion of organic nitrogen to a plant-available inorganic form and hence downregulating photosynthesis.
Enhanced R a in CTRL and in GEOENG with respect to BASE are caused by enhanced GPP, which is attributed to the CO 2 fertilization effect (see supplementary figure S7). Lower R a in GEOENG than in CTRL is primarily caused by the cooler climate in GEOENG. NPP, which is determined by GPP and R a (see equation (1)), increases in both CTRL and GEOENG compared to BASE (figures 3(g) and (h)). Stronger reductions in GPP than in R a result in lower NPP in GEOENG with respect to CTRL ( figure 3(i)). Like R a is sensitive to temperatures and GPP variations, R h is sensitive to changes in temperature and the carbon amounts in litter and soil. Higher rates of GPP in CTRL and GEOENG compared to BASE drive higher rates of R h because of larger litter inputs that also increased soil pools (see supplementary figure  S8). Therefore, higher R h is simulated in both CTRL and GEOENG compared to BASE in most regions (see supplementary figure S9). Nevertheless, regions with reduced precipitation and lower temperatures in GEOENG than in CTRL, such as the northern Amazon, undergo lower R h and evapotranspiration rates in GEOENG, which in turn retains more water in soil due to stomatal closure in plants (Swann et al 2016) (see supplementary figures S5 and S10).
The spatial patterns of NEP are determined by NPP and R h . Both CTRL and GEOENG demonstrate enhanced NEP compared to BASE with the exception of North America. Higher NPP induces higher NEP in CTRL and in GEOENG with respect to BASE; however, stronger increases of R h due to land use change and accelerated litter input than the increases of NPP in North America cause reduced NEP in both CTRL and GEOENG when compared to BASE (see supplementary figure S11). Similar to NPP, the cooler climate in GEOENG leads to lower NEP than that of CTRL. The variations of NEP and disturbance determine the perturbations of NBP, which represents the long-term and large-scale carbon uptake by ecosystems (see equation (1)). In GLENS outputs, disturbance is excluded in NEP but is included in NBP. Hence, differences in the spatial patterns between NBP and NEP illustrate the human and natural disturbances. In general, CTRL and GEOENG simulate enhanced NBP attributed to stronger GPP than in BASE (figures 3(j) and (k)); smaller NBP over North America in CTRL than in BASE is caused by increased carbon loss due to fires (see supplementary figures S12 and S13). Similarly, smaller burned area and lower carbon loss due to fires in GEOENG explain the enhanced NBP in North America, the boreal region of Eurasia, and the southeast coast of Australia with respect to CTRL (figure 3(l)).

Carbon sink strength and atmospheric CO 2 trajectory adjustments
The carbon sink strength (CSS) is determined by the accumulated NBP over a certain period of time. Table  1 lists the accumulated global carbon amount changes over various time periods for each carbon fluxes in Table 1. Accumulated global climate change over various time periods. Numbers represent the mean and cross-ensemble standard deviation of the first three ensemble members for the BASE, CTRL, and GEOENG (abbreviated as GEOG). PRECIP and ET are precipitation and evapotranspiration (unit: mm day −1 ), respectively, and the other variable annotations (unit: Pg C) are the same as equation (1) The trajectories of (a) accumulated global carbon sink strength changes (Pg C) due to geoengineering, (b) the atmospheric CO2 mole fraction (ppm) during 2010-2097 for BASE + CTRL (black line) and the adjusted atmospheric CO2 mole fraction due to terrestrial BGC feedbacks under geoengineering (blue line), (c) surface temperature responses (K) due to atmospheric CO2 adjustments, and (d) sulfur injection rates (Tg yr −1 ) in GLENS (red) and adjusted injection rates due to terrestrial BGC feedbacks (blue). equation (1). The terrestrial ecosystem has higher global total GPP primarily due to elevated CO 2 levels in both CTRL and GEOENG with respect to BASE. All the other terrestrial carbon fluxes except NBP in the terrestrial carbon budget (see equation (1)) are enhanced over time in both CTRL and GEOENG. Increases of GPP, NPP, NEP, R a , and R h in GEOENG are slower compared to CTRL because of lower temperatures and lower precipitation rates in GEOENG. Differences of accumulated NBP between GEOENG and CTRL show stronger CSS in GEOENG than in CTRL by 79 ± 6 Pg C at year 2097, increasing from 13 ± 2 Pg C during 2020-2039 to 30 ± 3 Pg C during 2078-2097. Such CSS enhancement is attributed primarily to reduced R a and R h , which leaves more carbon resident in land globally. Additionally, smaller carbon losses due to natural or anthropogenic disturbances (estimated by NEP−NBP) in GEOENG also indicate a larger land carbon reservoir.
For individual ecoregions, about 63% (125 ± 1 Pg C) of global NBP (198 ± 4 Pg C) in CTRL and 47% (130 ± 1 Pg C) of global NBP (277 ± 3 Pg C) in GEOENG during 2020-2097 are contained in evergreen broadleaf forests in tropical regions, which have the largest CSS among all ecoregions (see supplementary table S1). Nevertheless, the largest CSS differences (∆CSS) between GEOENG and BASE are found in croplands and mixed forests (see supplementary figure S14). The former are responsible for 27% (22 ± 2vPg C) and the latter for 25% (19 ± 0 Pg C) of the total ∆CSS (79 ± 6 Pg C) by 2097. Less heat stress due to lower temperatures in GEOENG than in CTRL is the main factor that increases crop yields (i.e. more carbon gains on land) in spite of lower precipitation (Proctor et al 2018). Increases in ∆CSS are also simulated in all other vegetated ecoregions during 2020-2097, primarily as a result of reduced respiration rates. Most ecoregions in both CTRL and GEOENG undergo increasing CSS from 2020-2039 to 2050-2069 but declined CSS from 2050-2069 to 2078-2097. Land use changes due to forest clearing and conversions to agriculture in the RCP8.5 emission scenario are likely the cause of these CSS changes (Hurtt et al 2011).
The changes in ∆CSS provide information to reconstruct the atmospheric CO 2 trajectory in GLENS experiments in which the same atmospheric CO 2 mole fractions based on the RCP8.5 emission scenario were used to drive the simulations. Even though carbon pools were simulated with an active terrestrial carbon cycle, terrestrial BGC feedbacks to the Earth system were not accounted for in the coupled modeling system. Such feedbacks can alter the carbon amount in the atmosphere, resulting in a different atmospheric CO 2 trajectory and hence climate change. Given the assumption that the land carbon fluxes in CTRL are consistent with the specified atmospheric CO 2 trajectory, we evaluate carbon gains on land as a result of aerosol geoengineering through comparing ∆CSS between GEOENG and CTRL (figure 4(a)). The global ∆CSS reaches 79 ± 6 Pg C by 2097. That is, terrestrial ecosystems globally store additional carbon in GEOENG compared to CTRL, amounting to 37 ± 3 ppm CO 2 -equivalent, reducing the atmospheric CO 2 mole fraction by as much as 4% in 2097 ( figure 4(b)). The reduced airborne fraction of anthropogenic carbon due to terrestrial BGC feedbacks is within the range reported in previous studies by the year 2100 despite different models and aerosol geoengineering strategies compared with GLENS (Keller et al 2014, Tjiputra et al 2016, Cao andJiang 2017).
By accounting for the terrestrial BGC feedback, lower atmospheric CO 2 levels would induce a cooler global surface temperature in GEOENG compared to the prescribed CO 2 trajectory, differing by 0.14 K in 2097 (figure 4(c)). Hence, in order to maintain the global surface temperature at 2020-levels, less sulfur aerosol would have needed to be injected (approximately 1.4 Tg yr −1 less at 2097) had we used prescribed emissions rather than prescribed concentrations in GLENS ( figure 4(d)). In addition, the global carbon budget assessment indicates that annual anthropogenic carbon emissions are stored 31% in land, 23% in the ocean, and 46% in the atmosphere (Friedlingstein et al 2019). Thus, the atmospheric CO 2 trajectory, altered slowly by additional carbon uptake on land, would effectively reduce ocean carbon uptake, thereby partially compensating for the reductions in atmospheric CO 2 mole fraction. In general, marine CSS is enhanced in a cooler climate that maintains higher CO 2 solubility and stronger overturning circulation in oceans; however, increasing anthropogenic ocean acidification would likely mediate marine CSS due to its detrimental effects on calcifying biota (Orr et al 2005, Fabry et al 2008, Doney et al 2009, Gattuso et al 2015. While fully coupled emissions-forced simulations with interactive terrestrial and marine biogeochemistry are required to quantify competing feedback effects, we expect that temperature-induced reduction of the solubility pump in CTRL simulations would be responsible for stronger marine CSS in GEOENG simulations; however, over decadal timescales (out to 2100), we expect a smaller ocean feedback compared to the land feedback in GEOENG, consistent with other studies (

Conclusion
This study characterizes responses and feedbacks of terrestrial ecosystems from an ensemble of coupled high-resolution Earth system model climate change simulations under a high greenhouse gas scenario with a geoengineering mitigation strategy that completely offsets temperature warming due to elevated anthropogenic CO 2 . From the GLENS simulation results, we find that terrestrial ecosystems could store an additional 79 Pg C on land globally by the end of the twenty-first century due to SO 2 injections in the lower stratosphere. Such carbon gains on land are mainly attributed to lower ecosystem respiration and diminished disturbance effects under geoengineering. As a result, less geoengineering effort would be needed to maintain the global surface temperature at 2020-levels in GLENS. Since GLENS simulated with an active terrestrial carbon cycle but without marine BGC feedbacks, we suggest that conducting fully coupled and emissions-forced ESM simulations with both marine and terrestrial BGC feedbacks enabled, as well as including stratospheric and tropospheric chemistry, will be necessary to reduce the uncertainty in quantifying aerosol geoengineering impacts on the Earth system.

Acknowledgments
This research was supported by the Reducing Uncertainties in Biogeochemical Interactions through Synthesis and Computation (RUBISCO) Science Focus Area (SFA), which is sponsored by the Regional and Global Model Analysis activity of the Earth & Environmental Systems Modeling Program in the Earth and Environmental Systems Sciences Division (EESSD) of the Office of Biological and Environmental Research (BER) in the US Department of Energy Office of Science. This research was also supported by the Oak Ridge National Laboratory Terrestrial Ecosystem Sciences (TES) SFA, sponsored by the Terrestrial Ecosystem Sciences program, which is also part of EESSD in BER in the US Department of Energy Office of Science. This research used resources of the Oak Ridge Leadership Computing Facility (OLCF) at Oak Ridge National Laboratory (ORNL), which is managed by UT-Battelle, LLC, for the US Department of Energy under Contract No. DE-AC05-00OR22725. Support for B.K. was provided in part by the National Science Foundation through agreement CBET-1931641, the Indiana University Environmental Resilience Institute, and the Prepared for Environmental Change Grand Challenge initiative. The Pacific Northwest National Laboratory is operated for the US Department of Energy by Battelle Memorial Institute under contract DE-AC05-76RL01830. Support for L.X. was provided by the National Science Foundation grant AGS-1617844. The CESM project is supported primarily by the National Science Foundation (NSF). This material is based work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the NSF under Cooperative Agreement No. 1852977. Computing and data storage resources, including the Cheyenne supercomputer (doi:10.5065/D6RX99HX), provided by the Computational and Information Systems Laboratory (CISL) at NCAR. All simulations were carried out on the Cheyenne high-performance computing platform https://www2.cisl.ucar.edu/usersupport/acknowledging-ncarcisl and are available to the community via the Earth System Grid.

Contributions
C-EY performed the analysis, generated graphical results, and prepared the manuscript in consultation with FH, JF, LX, and DR. FH provided impulse response function calculations. ST, DM, JR, MM, BK were the core team conducting the model simulations in GLENS. All authors contributed to group discussions and manuscript revisions.