Effects of 21st‐century climate, land use, and disturbances on ecosystem carbon balance in California

Abstract Terrestrial ecosystems are an important sink for atmospheric carbon dioxide (CO2), sequestering ~30% of annual anthropogenic emissions and slowing the rise of atmospheric CO2. However, the future direction and magnitude of the land sink is highly uncertain. We examined how historical and projected changes in climate, land use, and ecosystem disturbances affect the carbon balance of terrestrial ecosystems in California over the period 2001–2100. We modeled 32 unique scenarios, spanning 4 land use and 2 radiative forcing scenarios as simulated by four global climate models. Between 2001 and 2015, carbon storage in California's terrestrial ecosystems declined by −188.4 Tg C, with a mean annual flux ranging from a source of −89.8 Tg C/year to a sink of 60.1 Tg C/year. The large variability in the magnitude of the state's carbon source/sink was primarily attributable to interannual variability in weather and climate, which affected the rate of carbon uptake in vegetation and the rate of ecosystem respiration. Under nearly all future scenarios, carbon storage in terrestrial ecosystems was projected to decline, with an average loss of −9.4% (−432.3 Tg C) by the year 2100 from current stocks. However, uncertainty in the magnitude of carbon loss was high, with individual scenario projections ranging from −916.2 to 121.2 Tg C and was largely driven by differences in future climate conditions projected by climate models. Moving from a high to a low radiative forcing scenario reduced net ecosystem carbon loss by 21% and when combined with reductions in land‐use change (i.e., moving from a high to a low land‐use scenario), net carbon losses were reduced by 55% on average. However, reconciling large uncertainties associated with the effect of increasing atmospheric CO2 is needed to better constrain models used to establish baseline conditions from which ecosystem‐based climate mitigation strategies can be evaluated.


| INTRODUC TI ON
Changes in land use have been a primary factor in the global rise of atmospheric carbon dioxide (CO 2 ) and a major driver of global climate change (Le Quéré et al., 2017). Since 1850, landuse change (LUC) has added nearly half as much carbon to the atmosphere as fossil fuel emissions and has exerted a dominant influence on the storage of carbon in terrestrial ecosystems (Houghton & Nassikas, 2017;Le Quéré et al., 2017). Changes in land use have primarily resulted in the conversion of natural ecosystems to produce food, fiber, and for establishment of settlements (Foley et al., 2005). Combined, land use and management have the potential to undermine the ability of ecosystems to produce a wide range of services (Foley et al., 2005), including the storage and sequestration of carbon to mitigate climate change (Houghton & Nassikas, 2017). Meanwhile, climate change affects ecosystem carbon balance by changing the rate of carbon uptake in vegetation (Ballantyne, Alden, Miller, Tans, & White, 2012;Ballantyne et al., 2017) and the decay and decomposition of dead organic matter (DOM) and soils (Melillo et al., 2017;Pries, Castanha, Porras, & Torn, 2017), and remains the subject of intense study (Fahey, Doherty, Hibbard, Romanou, & Taylor, 2017;Domke et al., 2018). Furthermore, studies indicate that climate change is increasing the frequency and magnitude of extreme events (Mann et al., 2017) which can alter ecosystem carbon balance by increasing gaseous emissions and through the transfer of carbon from live to DOM pools (Kurz et al., 2008). The combined effects of climate and land change can result in either positive (i.e., net emissions) or negative (i.e., net sequestration) feedbacks from the biosphere on the climate system (USGCRP, 2017). The direction and magnitude of these feedbacks will either hinder or facilitate the achievement of local to global-scale greenhouse gas reduction targets.
Accounting for and minimizing anthropogenically linked ecosystem carbon emissions and/or maximizing the biosphere carbon sink is important for meeting the 2°C target of the Paris Agreement (Fargione et al., 2018). Despite the historically negative role in the carbon cycle, land use and management is increasingly being evaluated instead as a tool for climate mitigation (Cameron, Marvin, Remucal, & Passero, 2017;Fargione et al., 2018;Griscom et al., 2017;Houghton, Byers, & Nassikas, 2015). To evaluate the effectiveness of land-based mitigation efforts, regional baseline estimates of carbon stocks and fluxes are needed. At management scales, spatially explicit projections spanning long temporal periods (i.e., 50-100 years) need to consider the interactive effects of changes in land use/land cover (LULC) and climate. Such projections provide a reference point against which the effects of mitigation actions, applied to different locations at different times, can be evaluated. With the rise in subnational emission reduction targets (Galarraga, de Murieta, & França, 2017) jurisdictions must work not only to mitigate emissions from both the energy and land sectors but also to understand future climate-biosphere feedbacks and their effect on ecosystem carbon balance.
The State of California exemplifies many of the challenges associated with projecting integrated effects of climatic and LUC on ecosystem carbon balance. California is a large and ecologically diverse region overlain by a complex land-use mosaic (Figure 1).
California ranks first among the United States in population, economic activity, and agricultural production value, while at the same time maintains nearly half of its land area in some form of protected status. California is a leader in subnational climate action with one of the only economy-wide regulatory climate targets in the world (U.S. Climate Alliance, 2018). Recognizing the importance of ecosystems for emissions mitigation, the state has started to integrate ecosystem management as part of its emissions reduction strategy (California Air Resources Board, 2017 (Swain, Langenbrunner, Neelin, & Hall, 2018), with a recent severe drought killing more than 100 million trees (Stephens et al., 2018) and a high probability of a multidecadal drought occurring this century (Cook, Ault, & Smerdon, 2015). Anthropogenic climate change is increasing the frequency of climate extremes in California (Cvijanovic et al., 2017;Diffenbaugh, Swain, & Touma, 2015), leading to changes in natural disturbance regimes (Abatzoglou & Williams, 2016;Crockett & Westerling, 2018). However, the future effects of climate and land change on ecosystem carbon balance are difficult to predict, and pose additional challenges to meeting greenhouse gas reduction targets. Furthermore, no existing studies have attempted to quantitate the effects of major controlling processes such as land use, LUC, natural disturbance, and climate change, on the carbon balance of California ecosystems.
We consider 32 unique future scenarios exploring alternative assumptions about land use and LUC, reductions in global emissions, and future climate conditions. We used a fully coupled state-andtransition simulation model with carbon stocks and flows to estimate changes in ecosystem carbon balance for California's natural and agricultural lands on an annual time step for the period 2001-2100 (Daniel, Sleeter, Frid, & Fortin, 2018;Sleeter et al., 2018). Scenario projections explored low, medium, high, and business-as-usual (BAU) LUC pathways  and were combined with climate projections for two radiative forcing scenarios (i.e., representative concentration pathway [RCP] 4.5 and RCP 8.5). RCPs were simulated by four bias-corrected, statistically down-scaled global climate models (GCMs) (Pierce, Cayan, & Thrasher, 2014) chosen to represent a wide range of future climate pathways in California. Finally, we integrated projections of climate-driven drought-induced tree mortality and wildfire, overcoming a major limitation of previous studies (Allen, Breshears, & McDowell, 2015). The goals of this study were to (a) estimate changes in California ecosystem carbon balance and their uncertainties over the recent historical past and under a range of plausible future scenarios; (b) estimate how changes are distributed across different components (i.e., carbon pools), LULC classes (e.g., forests, grasslands), and regions; and (c) develop understanding of the relative impact of major controlling processes such as land use, disturbances, and climate change.

| MATERIAL S AND ME THODS
We used the land use and carbon scenario simulator (LUCAS), an empirical model of LUC coupled with a gain-loss model of ecosystem carbon dynamics (Selmants, Giardina, Jacobi, & Zhu, 2017;Sleeter et al., 2018; to project changes in ecosystem carbon balance for the state of California under a range of climate and land-use scenarios. The LUCAS model utilizes a fully coupled state-and-transition simulation model with stocks and flows (STSM-SF; Daniel et al., 2018) to estimate changes in ecosystem carbon pools resulting from changes in land use, land cover, and disturbances . The model estimates annual changes in carbon pools resulting from vegetation productivity, litterfall, mortality, decay/decomposition, emission, leaching, and harvest.
Carbon stocks and fluxes respond to changes in land use and land cover resulting from processes associated with urbanization, agricultural expansion and contraction, forest and agricultural harvest and management activities, and wildfire and drought-induced tree mortality. The model also considers the effects of short-and long-term climate variability on growth of live biomass and the turnover of DOM. Land-use transitions in this study are based on a set of projections developed for the state of California (Sleeter, Wilson, Wilson, Sharygin, & Sherba, 2017). Carbon stocks and fluxes are based on a national assessment of ecosystem carbon balance  modified for the State of California based on methods described in Daniel et al. (2018).

| Study area
The spatial extent of this study was the state of California covering 423,812 km 2 ( Figure 1). The State of California was subdivided into a regular grid of 1 km × 1 km simulation cells. The state type of each simulation cell was based on combinations of 12 ecological regions, 58 counties (administrative units), and 12 discrete LULC classes, including water, wetlands, snow/ice, barren, forest, grassland, shrubland, annual cropland, perennial cropland, development, and transportation/roads. For the forest, grassland, shrubland, and cropland classes we also tracked the age and time-since-transition (TST) of each cell as additional state variables. A general description of the structure of STSMs can be found in Daniel, Frid, Sleeter, and Fortin (2016).

| States and transitions
The STSM divides a landscape into a grid of simulation cells and then simulates the state type of each cell forward in time, as a discretetime stochastic process, in response to any number of possible transitions (Daniel et al., 2016). Transitions between state types were defined to represent the processes associated with urbanization, agricultural expansion and contraction, agricultural harvest and management, forest harvest, wildfire, and drought mortality ( Figure 2). In total, 39 unique transition pathways were defined for this study and were based on the model described in Sleeter, Wilson, et al. (2017) and . The order in which transitions were applied was randomized in each time step. Several important modifications to the model were made for this study and are described below. For a thorough description of the structure of the STSM developed for California, see Sleeter, Wilson, et al. (2017); and Wilson, Sleeter, and Cameron (2016).
The base land change model included a BAU and three population-based scenarios which explored alternative urbanization pathways (Sleeter, Wilson, et al., 2017;. The BAU scenario was based entirely on extending historical rates of LULC change into the future and represents a future scenario pathway describing a continuation of recent historical trends. The three population-based scenarios were originally designed to explore the effect of different rates of population growth on urban land expansion. For this study, we further modified the population-based scenarios to also include varying rates of other anthropogenic land uses, such as agricultural expansion and contraction and forest harvest.
The 'high' land-use scenario utilized a high population projection to characterize urbanization, and was further modified to explore a future with high rates of other anthropogenic land uses. Under this scenario, transition rates associated with agricultural expansion and contraction were derived by sampling, with replacement, from the period 1997-2002, which was a period of sustained growth in the areal extent of California agricultural lands (Sleeter, Wilson, et al., 2017;. Similarly, for transitions associated with forest harvest, we sampled, with replacement, from the 2002-2009 period. A similar approach was used to construct the 'low' land-use scenario, where the rates of agricultural land use were sampled from the 1993-1996 period and transitions associated with forest harvest were sampled from the 2010-2014 period. These rates were combined with a low population projection and corresponding low rate of urbanization. For the 'BAU' and 'medium' scenarios, we sampled, with replacement, from the full time series of historical land-use data for each transition type. As such, the only difference between the 'BAU' and 'medium' land-use scenarios was the rate of urbanization. The 'BAU' scenario projected urbanization using the full time series of historical urbanization data while the 'medium' land-use scenario used a medium population projection as described in Sleeter, Wilson, et al. (2017) and .
Forest harvest activities in California are dominated by softwoods, which make up approximately 58% of the total forest area of the state and are dominated by a mixed conifer group (includes mixed stands of Douglas fir, ponderosa pine, sugar pine, Jeffrey pine, white and red fir, incense cedar, and other true fir species), ponderosa pine, and other western softwoods (Christensen, Waddell, Stanton, & Kuegler, 2016). Hardwoods make up approximately 40% of the state's forest lands and are dominated by western oak species; however, they make up less than 1% of the total statewide harvest (McIver et al., 2015). Conversely, the combined harvest of Douglas fir, true firs, ponderosa pine, redwood, and sugar pine accounts for nearly 94% of all harvest in California (McIver et al., 2015). In the LUCAS model, forest harvest was characterized as either stand replacing (clear-cut) or partial harvest (selection) events.
Annual historical harvest rates were derived from annual historical  For selection harvest, we assumed 20% of live biomass pool was transferred to the harvested wood products (HWP) pool.
Wildfire and drought-induced tree mortality were estimated based on submodels which were integrated into the LUCAS frame- work. An exogenous statistical fire model was used to derive burn area projections for each climate scenario (Sleeter, Wilson, et al., 2017;. For the projected period, burn area was estimated for each GCM and RCP based on a statistical model of wildfire which considered the effects of climate, vegetation, population density, and fire history (Westerling, 2018 We fit a binomial generalized linear model for each of the three mortality classes, using SPEI as a single predictor in the model. We used these models to estimate future drought-induced mortality for each GCM and RCP scenario on an annual time step. We estimated the annual mortality area for each ecoregion from the model outputs and sampled, with replacement, from a Gaussian distribution created from these annual ecoregional means and an assumed 50% standard deviation. Annual relative probability maps derived from the spatial predictions were used within LUCAS to constrain the pattern of disturbance. Additional details describing the wildfire and drought-induced tree mortality submodels can be found in the Data S1: Methods.

| Carbon stocks and flows
With the STSM-SF method, the fate of continuous carbon stocks can be simulated for each simulation cell, based on a suite of continuous flows (i.e., carbon fluxes) specifying the rates at which these stock levels change over time . Carbon stocks considered in this study included a single live pool (living biomass), three DOM pools (standing deadwood, down deadwood, and litter), and a soil pool (SOC, soil organic carbon). In addition, we tracked carbon in three product pools, including HWP and carbon stored in agricultural products (grain and straw pools). Pools representing carbon stored in the atmosphere and aquatic pools were tracked to ensure the mass balance of carbon was maintained.
For live, DOM, and soil pools, initial carbon stocks were estimated for each land cover class (forest, grasslands, shrublands, and annual and perennial agriculture) and ecoregion based on a regionally calibrated dynamic global vegetation model (DGVM; Daniel et al., 2018;Selmants et al., 2017;Sleeter, Liu, Daniel, Frid, & Zhu, 2015;Sleeter et al., 2018) and a remote sensing-derived map of forest age (see Data S1: Methods).
In addition, LUCAS estimates change in carbon stocks and fluxes re- Ecoregion and land cover-specific turnover rates were then used to represent the flow of carbon from live to DOM pools resulting from litterfall and mortality based on an analysis of output from a DGVM Sleeter et al., 2018). Rates of decay and decomposition of DOM, as well as the lateral flux from of carbon from terrestrial to aquatic systems (i.e., leaching), were specified for each ecoregion ). An annual multiplier was applied to DOM turnover rates based on a Q10 function (Kurz et al., 2009). We

| Initial conditions
The LUCAS model was initialized with spatially explicit maps representing spatial strata, LULC type, and age ( Figure S1). All maps were projected to an Albers Equal Area Conical Projection using the NAD83 Datum. Each cell had a spatial resolution of 1 km × 1 km.  Wilson, et al. (2017) and  with modifications for perennial croplands described in Data S1: Methods.
Initial carbon stocks for each simulation cell were estimated based on the ecoregion, LULC class, and age of each cell and an age-to-carbon look-up table derived from the output of the DGVM Sleeter et al., 2018; Figure S2). Initial stock estimates were then scaled using a stationary spatial growth multiplier to reflect within-ecoregion heterogeneity of carbon pools. The spatial multiplier was estimated using 30 year climate normal's and the empirical NPP model, where values for each cell reflected the NPP anomaly relative to each cell's ecoregion and LULC class type. The soil carbon pool was estimated from a soil carbon stock spatial layer using soil property prediction maps (see Data S1: Methods). We used this soil carbon map to scale the regional output from the DGVM following methods developed and described in Daniel et al. (2018).

| Effect of CO 2 fertilization
Increasing atmospheric CO 2 (C a ) concentration has a direct, positive effect on carbon assimilation by plants through photosynthesis (Franks et al., 2013) where NPP t,s CFE is the projected NPP including CFE for year t (a year between 2001 and 2100) and scenario s (either RCP 4.5 or RCP 8.5), NPP t,s NOCFE is the NPP projected by the model in the absence of CFE for year t and scenario s, ΔC t,s a is the difference in C a between year t and the base year of 2001 for scenario s, and β is a percentage increase in NPP for a 100 ppm increase in C a .  (Table 1). During this period, the cumulative change in total ecosystem carbon was −163.1 Tg C, equivalent to approximately 34% of the state's total greenhouse gas emissions over  Table 1 shows the estimated carbon stocks and fluxes for each time step over the historical period.

| Contemporary change in carbon stocks
TA B L E 1 Annual carbon stocks and fluxes for California for the period 2001-2015. The transfer column represents the cumulative carbon transferred from ecosystem classes to nonecosystem classes (e.g., development). The LULCC column represents carbon losses due to land use, LUC, and disturbances

| Projected change in ecosystem carbon balance
This study projected changes in ecosystem carbon balance for 32 unique future pathways representing all possible combinations of two radiative forcing trajectories (RCP 4.5 and RCP 8.5) as simulated by four global climate models (CanESM2, CNRM-CM5, HadGEM2-ES, and MIROC5), and four land-use scenarios (BAU, high, medium, low;Sleeter, Wilson, et al., 2017). Each scenario was replicated for 100 Monte Carlo simulations. Figure 4 shows the mean projected total ecosystem carbon (Figure 4a), live, DOM, and soil carbon stocks (Figure 4b), and NECB ( Figure 4c)

| Effects of different scenario pathways
Climate models were the primary source of uncertainty in projected ecosystem carbon balance (Table 3; Figure 5). When averaged across all scenario simulations, mean annualized net carbon flux ranged from −8.4 to −0. Increases in precipitation under the CanESM2 RCP 8.5 scenario resulted in large increases in NPP which led to increased rates of sequestration not realized under the RCP 4.5 future (Table 3). These results suggest the overall benefits of global-scale reductions in radiative forcing and reductions in atmospheric greenhouse gasses may not be realized at local-to-regional scales and may result in positive feedbacks on the climate system through reductions in ecosystem sequestration potential.
The four land-use scenarios used in this study were developed to explore alternative future pathways based on low, medium, and high population and land-use intensity assumptions, along with a scenario which explored a continuation of current trends in land use and LUC (e.g., BAU; Sleeter, Wilson, et al., 2017;. Projected changes in land use are shown in Figure S3.
Urbanization was projected to occur primarily on agricultural lands, grasslands, and shrublands while changes in agriculture were primarily located in grassland and shrubland ecosystems. The 'high' scenario was the only scenario where agricultural lands were projected to increase in overall area ( Figure S4). Under all scenarios, forest, grassland, and shrubland were projected to decline from

| Effect of CO 2 fertilization
The effect of C a enrichment on ecosystem productivity has a major impact on uncertainty in the estimation of regional scale carbon balance (Figure 7). The inclusion of a CFE effect with β values ranging from 0.027 to 0.136 resulted in increases in carbon sequestration by the end-of-century ranging from 1.7-12.4 Tg C/year (assuming the CFE effect was limited to C a of 600 ppm). Assuming no CFE limitation based on C a , NECB was projected to increase between 5.7 and 24.3 Tg C/year under the RCP 8.5 scenario.
Compared to the reference scenario (β = 0), total carbon stored in California ecosystems was projected to increase by 12.5% under However, uncertainties resulting from the effect of increases in atmospheric CO 2 are large, and must be reconciled in order to better constrain model projections.

| Implications of different climate and landuse pathways
By comparing various scenario combinations, we can determine the relative impact of both global adherence to a lower emission trajectory and jurisdictional actions to reduce emissions from land use and LUC (Figure 6). On average, the cumulative effect of a

| Driving forces of carbon change
In this study, we estimated the direct effects of climate change on ecosystem carbon balance through both changes in vegetation productivity and the decay and decomposition of DOM. Additionally, projected changes in climate were used to estimate changes in the magnitude and frequency of natural disturbance events including wildfire and drought-induced tree mortality. The climate models considered in this study provided a wide range of future conditions, with mean annual temperature increases ranging from 1.7 to 4.5 C by the end-of-century ( Figure S5). The effect of increased warming resulted in large declines in soil carbon pools across almost all scenarios ( Figure 4) GCMs show an upward trend in projected wildfire emissions, the much higher emissions from CanESM2 were due to a steep increase in wildfire area during 2080-2100 ( Figure S6). High-and mediumseverity wildfire also led to an additional pulse of carbon from live to DOM pools due to wildfire mortality. The effect of climate change on drought-induced tree mortality was only evident in the high-severity mortality class under RCP 8.5. High-severity mortality led to a projected cumulative average transfer of 613 Tg C from the live to the DOM pools. In all scenarios, extreme episodic mortality events were driving carbon losses rather than consistent low-level background events ( Figure S6). This is characteristic of severe drought periods driven by a combination of low precipitation and high temperatures, or high temperatures alone. Under RCP 4.5, high-severity mortality led to average cumulative losses of 563 Tg C from the live biomass pool, but no GCM differed more than 6% from this mean value.

| Comparison to other studies in California
The comprehensive nature of this study, which used a gain-loss method to estimate carbon stocks and fluxes for California's forest, grassland, shrubland, and agricultural ecosystems over both recent historical conditions as well as 32 alternative future pathways, represents a unique contribution to understanding the interactions and major controlling processes of ecosystem carbon balance. As such, we can compare our estimates of carbon stocks and fluxes with other recent studies which have been more limited in scope, often covering only a subset of the variables included in this study.
A recent study (Gonzalez, Battles, Collins, Robards, & Saah, 2015) estimated that in 2001, California ecosystems stored 920 Tg C in live aboveground vegetation, of which 830 Tg C was stored on 125,000 km 2 of land classified as forest. Our study estimated California forests stored 1638.8 Tg C in total live biomass carbon (sum of above-and belowground stocks). While we did not explicitly model above and below ground carbon stocks as separate pools, we can infer the aboveground portion by applying a standard root:shoot relationship of 0.27 (Mokany, Raison, & Prokushkin, 2006).
This provides an estimate of forest live aboveground carbon storage of 1,196.3 Tg C in 2001. Comparisons are further complicated due to differences in land classification, particularly differences in the extent and amount of forest area. We estimated the 2001 forested land area in California was ~156,000 km 2 , which was 20% higher than the estimate by Gonzalez et al. (2015). Normalizing for differences in the extent of forest area provides an estimate of 957 Tg C stored in forest aboveground live stocks, an amount more comparable to that provided by Gonzalez et al. (2015) (830 Tg C). For forests, the stock change study (Gonzalez et al., 2015) estimated a decline of −50 Tg C (−5.6 Tg C/year) between 2001 and 2010 which was approximately two times larger than our estimate of −25.6 Tg C (−2.8 Tg C/year) over the same period. However, the stock change approach did not include carbon accumulation (e.g., growth) on forest cells which did not change from one ordinal forest class (e.g., forest type, cover, and height) to another (Gonzalez et al., 2015), whereas this study's gainloss approach estimated growth on an annual basis for all simulation cells. A comparison of this study to other published research estimating ecosystem carbon stocks in California is shown in Table S1.
Wildfire is a major source of carbon emissions to the atmosphere in  (Reinhardt, Keane, & Brown, 1997) with fire footprints from CalFIRE (2016). Estimates from our study compare well with those from ARB over the overlapping historical period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) and are shown in Figure S9. We estimate that historical fire emissions averaged 14.3 (σ = 12.6) CO 2 per year compared to an average historical rate of 15.5 (σ = 11.2) CO 2 per year based on the ARB methodology. Annual estimates and variability also show a strong agreement with estimates produced using the ARB modeling approach with an r 2 value of 0.92 based on a comparison of annual emission estimates.

| Limitations of study
This study was limited to just 4 of the 32 LOCA down-scaled global climate models from the CMIP5 archive. While these four models were chosen by a working group for the California 4th Climate Assessment as representing a broad range of future California climates, our results would likely be different if all 32 models were assessed in this framework. More importantly, our current framework did not include variability in key parameters that would likely increase the uncertainty of the results. This includes the uncertainty associated with parameters related to climate (i.e., effect of precipitation and temperature on NPP and DOM turnover rates), carbon (base stocks and flows), and disturbance (extent of wildfire and mortality). Our approach did not account for changes in vegetation type which may result from the coupled effects of climate change and high-severity fire. Specifically, there is growing concern that in some regions, the ability of forests to recover after large high-severity fires is reduced in response to climate change (Thorne et al., 2017).
Our model did not estimate these effects and may result in an overestimation of carbon storage given the assumption that forests will always recover after disturbance.
Future warming, and its effect on DOM turnover rates, was represented using climate model temperature projections and a Q10 function generally consistent with those used in the Carbon While rising atmospheric CO 2 undoubtedly plays a role in vegetation productivity, the magnitude and duration of this effect and how it manifests in specific ecosystems is highly uncertain (Smith et al., 2016), especially under long-term scenario projections (Arora et al., 2013;Friedlingstein et al., 2014;Piao et al., 2013;Sitch et al., 2008). Early results from FACE studies appeared to support this strong positive CFE on NPP at the ecosystem level across a range of forest types (Norby et al., 2005). However, there is also substantial evidence that an array of feedback responses and constraints can dampen or eliminate CO 2induced increases in plant carbon assimilation and growth (Franks et al., 2013;Norby et al., 2010;Smith & Dukes, 2013). Longer term results from FACE studies in both forests and grasslands have demonstrated that the initial CFE-induced growth stimulation declined over time as nitrogen (N) became more limiting (Norby et al., 2010;Reich et al., 2006). The magnitude of this progressive N limitation on CFE is highly site-specific, and sometimes does not appear to dampen CFE at all because of high site N status or positive feedbacks to N cycling (Smith & Dukes, 2013 (Smith et al., 2016).
The β values included in the sensitivity analysis generally cover the range of values observed at forested FACE sites (Norby et al., 2005;Smith et al., 2016;Zaehle et al., 2014). These estimates tend to be higher than β values derived from a remote sensing-based approach (Smith et al., 2016), which have been criticized as likely and flux (NECB) is likely to be much larger than the estimates shown in Figure 7. The large degree of uncertainty resulting from the introduction of a range of CFE's in this study underscores the need for more research on the role of increasing atmospheric CO 2 on terrestrial carbon cycling, especially under long-term projections, as studies have indicated that the role of increasing CO 2 on ecosystem carbon balance can be as much as four times larger than the effect of climate change (Arora et al., 2013;Gregory, Jones, Cadule, & Friedlingstein, 2009). In the context of climate mitigation planning, assuming no CFE is the most conservative approach and that any realized CFE in the future will likely only further increase carbon sequestration in ecosystems.
Similar to most other carbon cycle models, our model might overestimate recovery rates of carbon sink strength after prolonged drought (Anderegg et al., 2015). Between 2013 and 2015, California entered a period of sustained and exceptional drought which affected most areas and ecosystems of the State. Evidence suggests the extreme drought conditions experienced in recent years were without precedent over instrumented and millennial records (Robeson, 2015), and while primarily the result of natural climate variability (Berg & Hall, 2015;Williams et al., 2015), anthropogenic warming was an import contributing factor (Diffenbaugh et al., 2015;Shukla, Safeeq, AghaKouchak, Guan, & Funk, 2015;Williams et al., 2015).  , 1976-1977, 1987-1992).
While the uncertainty associated with projecting extreme events is large (Swain et al., 2014), substantial evidence suggests anthropogenic climate change will increase the probability of extreme drought conditions in California Diffenbaugh et al., 2015;Williams et al., 2015) which may have significant impacts on ecosystem carbon balance.

| CON CLUDING REMARK S
Changes in ecosystem carbon balance are an important driver of global climate change. However, there is a great deal of uncertainty in the direction and magnitude of the carbon source or sink.
Results of this study show that uncertainty primarily resulted from future climate conditions driven by different climate models and are consistent with other studies in this regard (Zaehle et al., 2007). However, these uncertainties are relatively small compared to the uncertainties associated with increasing CO 2 and its effect on carbon storage and flux. While NECB was projected to increase relative to recent historical rates, carbon stored in California ecosystems was projected to decline in 30 of the 32 future scenarios considered in this study. Reducing global greenhouse gas emissions did not always result in increased carbon storage in California ecosystems, suggesting there may be unanticipated feedbacks to the climate system resulting from ecosystem carbon flux. Conversely, reducing LUC was a reliable and consistent approach to increasing carbon sequestration, regardless of future climate conditions.
These findings are useful for establishing a set of baseline projections from which additional ecosystem-based climate mitigation strategies can be evaluated (Cameron et al., 2017;Fargione et al., 2018). Such studies are needed to understand the role of ecosystems in regulating greenhouse gas emissions and for achieving local to global-scale climate mitigation goals.

ACK N OWLED G EM ENTS
The authors thank Zhiliang Zhu, Brad Reed, and Deb Willard of the U.S. Geological Survey for their ongoing support for this research. The authors also thank David Stoms of the California Energy Commission for his early guidance and support of land change scenario projections in California.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

CO D E A N D DATA AVA I L A B I LIT Y
The model, source code, and data required to replicate this study, as well as the output data supporting the conclusions of the study are available through the USGS ScienceBase repository.
The LUCAS model runs within the Syncro-Sim software application. All simulations were run using Syncro-Sim software version 2.0.18, under the Mono framework for Linux. The STSM and SF modules used in this study were version 3.1.18. We ran the simulations on the Comet system at the San Diego Supercomputing Center under the NSF Extreme Science and Engineering Discovery Environment (XSEDE) program through allocation TG-DEB17001767.
The following steps are required to run the model used in this analysis: • Download and install the latest Windows or LINUX version of the Syncro-Sim software, available at http://www.apexr ms.com.
• Download and unzip the model files, including spatial input files and a '.ssim' (SQLite) database from the online repository (to be filled in upon publication).
• Use the Syncro-Sim software to open the 'California Carbon Model.ssim' file and select a scenario to run.
• Alternatively, the model can be built from scratch using the R programming language with the rsyncrosim package installed. To follow this approach, download the 'California Carbon Model R Code.zip' data package and run the necessary R scripts.
All output results from the 32 scenarios described in this report are available from the ScienceBase online repository. Each scenario includes an SQLite database (SyncroSim '.ssim' file) and a compressed folder with all spatial output. Tabular