We classified climatic drylands based on the aridity index, using a threshold of <0.65 for the ratio of mean annual precipitation (P) to mean annual potential evapotranspiration (PET) (Yao et al., 2020) calculated using the 2.5 arc minutes 1981–2010 TerraClimate mean gridded surface climatology (Abatzoglou et al., 2018). The fine grain of the TerraClimate product corresponded more closely to physical reality than coarser resolution products (Supplementary Table S1 for details of these gridded products). We excluded drylands >55° north and south of the Equator to omit “cold” permafrost drylands, resulting in a global dryland area of 59.1 × 10⁶ km². While LSMs do not always simulate drylands exactly coincident with this climatic mask, this masking approach is standard (e.g., Ahlström et al., 2015; Yao et al., 2020; Gonsamo et al., 2021).

We used the PML-v2 (v016) GPP product derived using a light use efficiency model based on MODIS observations (Zhang et al., 2019a). This product has a spatial resolution of 500 m, revisit frequency of 8 days, and was available from 26/03/2000 to 26/12/2020. Importantly, the PML-v2 product couples evapotranspiration and GPP resulting in a more robust estimation of GPP while partly accounting for water use efficiency, and explicitly accounts for the influence of changing CO₂ concentration on carbon assimilation via a simplified photosynthesis model (Gan et al., 2018). PML-v2 has been evaluated at 95 eddy covariance flux tower sites, 40 of which lie within our delineated drylands (Zhang et al., 2019a). Over the 95 sites, PML-v2 outperformed other GPP products available at the time for 8-day and site mean GPP, with smaller bias, lower RMSE, and higher R² compared to the MOD17A2H, FluxCom GPP, and VPM GPP products; (Zhang et al., 2019a), while performance regarding annual anomalies was comparable. We calculated annual integrated productivity for spatial comparison with the LSM predictions and temporal analysis.

We inferred AGC in biomass from L-VOD retrieved from the Soil Moisture and Ocean Salinity (SMOS) mission L-band observations from 2011 to 2018 (SMOS-IC V2; Wigneron et al., 2021; Figure 1). Soil moisture and L-VOD are derived from a two-parameter inversion of the L-MEB model (L-band microwave emission of the biosphere) from the multi-angular and dual-polarized SMOS observations (Wigneron et al., 2021). The SMOS-IC products are mostly independent of other Earth observation datasets or simulations from atmospheric models (Wigneron et al., 2021). The L-VOD product consists of ascending and descending orbit datasets with a maximum 3-day revisit time. L-VOD is sensitive to radio frequency interference (RFI) over some specific geographic regions (particularly central Asia and southern Europe) and noise when retrieved over areas of variable topography or frozen ground (Fernandez-Moran et al., 2017). L-VOD data require filtering to retrieve reproducible signals, so we filtered observations by quality flags to exclude pixels containing more than 10% water, ice, or urban land cover (based on the MODIS MCD12Q1 product and IGBP classification scheme), extreme topography (Mialon et al., 2008) and frozen ground (ECMWF soil temperature <273 K) (Fernandez-Moran et al., 2017). L-VOD retrievals above a noise threshold between measured and L-MEB modelled brightness temperature values (8 K) were discarded to remove strong RFI (Fernandez-Moran et al., 2017). We selected four regions with sufficient data coverage (North America, South America, Australia, and Africa) and total dryland area 37.9 × 10⁶ km², as RFI precluded reliable retrieval of L-VOD in other dryland regions. Filtered observations (mean of 102 per pixel annually) were aggregated to annual median ascending and descending L-VOD. Where fewer than 20 observations per year were available, the pixel was excluded (Brandt et al., 2018). Large differences between annual average ascending and descending L-VOD can indicate remaining RFI, therefore pixels where this difference was greater than 0.05 L-VOD were excluded. This filtering also excluded a part of North America due to an L-VOD anomaly in 2011. L-VOD retrievals are sensitive to both vegetation biomass and water stress (Konings et al., 2019). However, we expect trends in annual median L-VOD to be largely unaffected by inter-annual variability in vegetation water content. We assume that any trend in vegetation water content would correspond to a trend in vegetation biomass in these predominantly water-limited ecosystems (Abel et al., 2021; Frappart et al., 2020). We assumed a linear relationship between plant water content and biomass in these generally sparsely vegetated drylands, which was supported by the findings of (Brandt et al., 2019; Fan et al., 2019; Konings et al., 2017; Liu et al., 2015; Mialon et al., 2020; Rodríguez-Fernández et al., 2018; Tian et al., 2016). Annual average carbon density (Mg C ha⁻¹) was estimated from ascending L-VOD, which is acquired at dawn when plant water storage usually peaks (Tian et al., 2018). Equation 1 was used to estimate AGC which was derived from an OLS regression of L-VOD against a global biomass map (Santoro et al., 2018), converted to AGC following the 47% biomass:carbon density ratio (Paustian et al., 2006; Chave et al., 2019). We developed a drylands-specific biomass transfer function because previous linear models (e.g., Brandt et al., 2018) calibrated against data including high-biomass ecosystems like tropical forests tend to overestimate dryland biomass (Supplementary Figure S3). The GlobBiomass map was selected as it showed a higher correlation with L-VOD for drylands than available alternatives (Supplementary Table S2). A G C = 52.48 × L V O D       [ M g   C   h a -1 ] (1)

To test the correspondence of state-of-the-art LSMs with retrieved biomass and productivity estimates, we used the simulation outputs from 12 LSMs from the TRENDY v8 project (Friedlingstein et al., 2019). This TRENDY ensemble of models includes CABLE-POP, CLASS-CTEM, CLM5.0, DLEM, ISBA-CTRIP, ISAM, JSBACH, JULES-ES-1.0, LPJ-GUESS, OCN, ORCHIDEE, and ORCHIDEE-CNP (Supplementary Table S3 for details). For these comparisons we focused on the most comprehensive simulations with time-varying CO_2, observed climate and land-use forcing (TRENDY “S3” simulations). Using the ensemble mean of these LSMs allows us to filter some of the inter-model variability to examine the overall model estimates. Six of these models include representation of fire processes [Supplementary Table S3 and (Friedlingstein et al., 2019) for details]. We annually aggregated GPP and mean annual biomass of the LSMs, from 01/01/2001 to 31/12/2018 and 01/01/2011 to 31/12/2018 respectively. We used AGC where this was quantified by the LSMs (CABLE-POP, CLASS-CTEM, CLM, ISAM). For the LSMs that did not explicitly partition carbon into above and below pools in TRENDYv8 outputs, we assumed 40% of simulated total biomass carbon was aboveground, using the average root to shoot ratios of grasslands, shrublands, savannas and woody savanna biomes that together dominate total dryland biomass (Liu et al., 2015; Qi et al., 2019). Though this ratio will vary regionally according to environmental factors and plant functional types, there is insufficient information available to enable accurate spatially explicit exploration.

To assess the implications of different model simulations over climate-relevant time scales, we evaluated changes in predicted soil C (cSoil, including necromass as litter and coarse woody debris where simulated), vegetation carbon (cVeg), and net ecosystem (cEco = cSoil + cVeg) stocks from 1901 to 2018. At the end of this period, we quantified how much each model diverged from the TRENDY mean in units of standard deviation of the ensemble mean (Goldstein et al., 2013). We excluded ISBA-CTRIP from this portion of the analysis because this model had issues with 1) insufficient soil carbon storage under trees and 2) high productivity of crops causing excessive soil inputs (Boysen et al., 2020) (subsequently resolved by the addition of a crop sub-model). To diagnose the causes of the differences in simulated carbon stocks over time and between models, we isolated the effects of CO₂, climate change, and land-use change by looking at the contribution associated with the TRENDY “S1” (including CO₂ forcing), “S2” (including CO₂ and climate change forcing), and “S3” simulations (Friedlingstein et al., 2019).

To understand the extent to which model-data differences in GPP and AGC related to fire, a process not represented in all LSMs, we investigated these spatially in relation to fire frequency. Using the most recent MCD64A1 version 6 burned area product (Giglio et al., 2015), we calculated the sum of the burned area between 2001 and 2018 per 1° grid cell and divided it by the length of the data record. This period was chosen to best encapsulate the influence of fires on GPP and AGC.

We tested the spatial correspondence between LSM simulations and biomass carbon retrieved from L-VOD or GPP retrieved from MODIS respectively. To maintain a consistent sample size across models, we used bilinear interpolation to resample the outputs from TRENDY models to a common spatial resolution of 1°. To appropriately handle boundary effects and partial coverage due to quality filtering of satellite retrieved L-VOD and GPP, we computed weights for the model pixels using the “exactextractr” package in R (Baston, 2020, v0.5.1). We computed weighted pairwise Pearson’s r (R package “weights”, Pasek, 2020, v1.0.1), and used weighted total least squares regression (R package “deming”, Therneau, 2018, v1.4) to fit linear models that account for uncertainty on both axes to quantify the agreement as the slope of the fitted model.

To test the temporal correspondence at annual resolution between biomass carbon inferred from L-VOD and simulated with LSMs and between GPP inferred from MODIS and simulated with LSMs, we summed the weighted dryland values at the native spatial resolution of each dataset using “exactextractr”, and for the biomass comparison considered only the areas with reliable L-VOD data for all years (Figure 1 and Supplementary Figure S2). We quantified bias as the mean error between the summed model predictions and the satellite retrievals. Both GPP and AGC measurements were highly sensitive to the approaches used to resample data. Sensitivity analysis revealed that the inclusion or exclusion of cells only partially within our spatial region of interest potentially introduces a 6-fold difference in the total productivity and biomass retrieved for drylands (Supplementary Figures S17–19). To minimise this issue when undertaking this analysis across datasets with different native resolutions, we used weighted extraction methods to account for the partial coverage of cells within a mask.

We used a Theil-Sen estimator to robustly fit linear models to the time series of productivity and AGC (Myers-Smith et al., 2020). To assess correspondence in inter-annual variability, we normalised each time series to its mean and detrended the series using the Theil-Sen model slope. The goodness-of-fit was then quantified as the mean absolute error between the series.

We found good spatial correlation between the mean GPP PML-v2 inferred from MODIS observations and GPP simulated by the TRENDY model ensemble (Figure 2; Table 1) over arid and semi-arid regions between 2001 and 2018. Correlation coefficients were moderately high across the 12 models, with Pearson’s r ranging from 0.71 to 0.90. The spatial patterns shown in Figure 2B are common across most LSMs (Supplementary Figure S6), with overestimation of productivity in the African Sahel and underestimation of productivity in southern African and South American drylands. The linear model slopes between the model and the satellite retrievals ranged from 0.56 to 1.18. The TRENDY mean exhibited strong correspondence with the satellite retrievals, with r of 0.9 and a slope of 0.88. We found no overarching relationship between fire frequency and bias in GPP although several models (JULES, CLASS-CTEM, ORCHIDEE-CNP, and OCN) did systematically overestimate GPP in more frequently burned areas (Supplementary Figure S11).

We found strong correspondence in inter-annual anomalies in GPP between MODIS-derived and all TRENDY models between 2001 and 2018. There was relatively little bias in mean GPP between the TRENDY mean and MODIS PML-v2 product (ME −1.5 Pg C yr⁻¹, 5.7%) (Figure 3C; Table 1). All 12 models and consequently the TRENDY-mean exhibited positive trends in productivity over time (11 of 12 trends were significant at α 0.05), while the retrieved GPP had no trend over the 18 years (Figure 3A, Supplementary Table S4). GPP trends differed spatially, increasing in some areas and decreasing in others (Supplementary Figures S5, S7).

We found often poor agreement in the spatial patterns of remotely-sensed and model-simulated biomass for the regions analysed (Africa, Australia, North and South America) (Figure 4, Table 1). Correlation coefficients were moderate across the 12 models, with Pearson’s r ranging from 0.52 to 0.75. The bias between the models and the satellite retrievals was often extreme, with most models exhibiting substantial bias and linear model slopes ranging between 0.2 and 5.2. The TRENDY mean exhibited slightly stronger correspondence with the satellite retrievals, with r of 0.75 and a slope of 1.13. For most models, we found little relation between observed burn frequency and model - data residuals, apart from JULES and CLASS-CTEM which overestimated biomass in more frequently burned regions (Supplementary Figure S12).

Detection of temporal change in biomass was hindered by the short (2011–2018) period included in the L-VOD product. There was generally poor agreement in the temporal patterns of simulated and remotely-sensed above ground carbon. There was a large range in bias between models and the satellite retrieval with MEs from −8.7 to 8.4 Pg C (−84.4 and 81.5%). Although inter-annual anomalies in AGC generally had the same sign across the L-VOD retrieval and TRENDY models (for the TRENDY mean the sign of the anomaly was consistent between TRENDY and the L-VOD in 88% of the time series), most LSMs had much smaller anomalies from the mean compared to the larger anomalies in the L-VOD-derived values (Figure 5). There were no meaningful trends in the time series of AGC; the few that were statistically significant had negligible effect sizes (Supplementary Table S4).

The land surface models predicted a net gain of ecosystem carbon in global drylands (mean 3.37 Pg C), due to an increase in soil carbon (mean 6.35 Pg C) that exceeds the losses in vegetation carbon (mean −2.98 Pg C) over the last 118 years (1901–2018) (Figures 6, 7, Supplementary Table S5). There was a wide range of carbon accumulation responses across the TRENDY models (Supplementary Table S5). The overall loss in simulated vegetation carbon was predominantly driven by the influence of land-use change, exceeding the increase simulated due to CO₂ fertilization (Supplementary Figure S13). Climate change had minimal effect on the simulated dryland vegetation carbon stocks. The overall increase in simulated soil carbon was predominantly driven by the increase in litter inputs into the soil associated with the increase in net primary productivity due to CO₂ fertilization, with most models predicting a small reduction due to climate change (Supplementary Figure S14). The influence of land-use change on soil carbon was more varied between model simulations, with the majority showing a reduction. The overall increase in simulated ecosystem carbon was predominantly driven by the influence of CO₂ fertilization, exceeding the losses predicted due to both land-use change and, to a lesser extent, climate change (Figure 6 and Supplementary Figure S15). Overall, dryland carbon stocks are simulated by models to have decreased by an average of 3.82 Pg C from 1901 to late-1960, before subsequently increasing by 7.19 Pg C up to 2018 as the CO₂ fertilisation effect increases simulated GPP in response to increasing atmospheric CO₂ concentrations.

Our results suggest that carbon uptake (GPP) by dryland ecosystems is broadly well represented by LSMs at the global scale. We found good agreement in the spatial distribution of GPP between simulated LSM estimates and MODIS based PML-v2 values in drylands for most models. Correlation coefficients were high and bias was low, especially for the TRENDY ensemble average. Our finding that spatial patterns of GPP were broadly in agreement was consistent with previous evaluations of other MODIS-derived GPP products against GPP observations from dryland eddy covariance flux towers in North America (Biederman et al., 2017). Geographic biases were consistent across most LSMs, with overestimation of productivity in the African Sahel and underestimation of productivity in drylands in southern Africa and South America.

We found very strong agreement in GPP interannual variability between LSM-simulated values and those modelled based on satellite observations of photosynthetic capacity (Figure 3B; Table 1), including the large 2 Pg C yr⁻¹ positive anomaly in 2010–11 caused by a severe La Niña event (Poulter et al., 2014; Haverd et al., 2016). This agreement lends some confidence that globally, LSMs capture the predominant dynamic responses controlling inputs of carbon into dryland ecosystems in response to annual-scale perturbations. However, the agreement found in IAV in GPP at the global scale contrasts with the disagreement in IAV at the site level where both MODIS and LSM simulated productivity have been found to underestimate IAV. Using eddy covariance observations from 25 dryland sites in North America, Biederman et al. (2017) found IAV in GPP was underestimated by MODIS-derived GPP. Similarly using observations from 12 eddy covariance sites in the southwestern U.S., MacBean et al. (2021) concluded LSMs underestimate IAV in dryland net ecosystem exchange. These findings underscore the ongoing need for further development and evaluations of remotely sensed GPP products in drylands (Smith et al., 2019).

All 12 LSMs consistently estimated a positive trend in dryland GPP between 2001 and 2018 (Figure 3; Table 1), mainly due to simulated CO₂ fertilization (Supplementary Figure S16). The CO₂ fertilization effect has been identified for drylands using satellite measurements (Donohue et al., 2013). Models predict a direct CO₂ fertilization effect (reduced photorespiration) and in water-limited systems stomatal closure for the same gain in CO₂ implies a longer growing season and higher annual productivity (Gonsamo et al., 2021), as expected from theory (Prentice et al., 2001). However, there was no trend in the satellite retrieved PML-v2 GPP product over this time. The discrepancy in GPP trends between models and PML-v2 was most pronounced in the Sahel region, where PML-v2 GPP decreased and LSM-simulated GPP increased (Supplementary Figure S7). While the PML-v2 product represents the effects of increasing CO₂ concentration on GPP (Zhang et al., 2019a), a comprehensive review of the effects of CO₂ enrichment on GPP by Walker et al. (2021) suggests that most remotely sensed GPP products are insufficiently sensitive to increasing CO₂. Satellite-retrieved GPP products are derived from light use efficiency models which have limitations (Smith et al., 2019; Baldocchi, 2020; Walker et al., 2021), one being the often very heterogeneous patterns of dryland vegetation and soil background signal giving rise to uncertainties, particularly in leaf area index estimations at the coarser spatial resolutions of the satellite pixels (Smith et al., 2019). Furthermore, the observations of evapotranspiration, energy, and carbon fluxes used for calibrating satellite retrieved GPP products and some of the LSMs are observed at relatively few locations globally, with short time-series, insufficient replication, and sampling biases that under-represent highly dynamic dryland ecosystems (Schimel et al., 2015; Hill et al., 2017; Jung et al., 2020). The undersampling of semi-arid settings in particular has been suggested to have a large impact on GPP upscaling errors (Jung et al., 2020).

Although we found reasonably good spatial agreement in AGC between L-VOD retrieved and average LSM (TRENDY-mean, Pearson’s r: 0.75), individual models exhibited poor agreement both with each other and with the satellite-derived estimates (Figure 4; Supplementary Figure S8; Table 1). We found large biases in total AGC between models (Figure 5A). To some extent, inter-LSM differences arise from different representations of land surface processes (Sitch et al., 2008), for example, their representation of fire and land-use change, and native model resolutions that influence PFT fractions within grid-cells along dryland fringes or vegetation transition zones (Supplementary Figures S17–19). These factors also influence simulated residence times (Friend et al., 2014) and ecosystem respiration, explaining why LSMs that show better correspondence with the PML-v2 product for GPP do not necessarily correspond better for AGC (Table 1).

To improve the evaluation of LSM vegetation carbon stock predictions in low biomass drylands, further efforts are needed to improve the accuracy and validation of satellite-derived AGC products. A recent review concluded existing aboveground biomass products are almost entirely inconsistent across dry forests, savannas, and grasslands (Zhang et al., 2019b). Uncertainties arise from insufficient data on non-forest vegetation for calibration and validation as well as the insensitivity of remotely sensed observations to low biomass ecosystems (Duncanson et al., 2019; Cunliffe et al., 2021). Furthermore, comparisons between LSM predictions and AGC maps are strongly influenced by root-shoot ratios. These ratios are poorly constrained by observations, particularly in drylands (Qi et al., 2019), and also vary under different environmental conditions in response to differences in atmospheric CO₂, aridity and grazing pressure, amongst other factors (Mokany et al., 2006; Yan et al., 2020; Yan et al., 2021). The impact of different root-shoot ratios also relates to the functional role of roots in LSMs. For many models, root biomass does not impact function such as explicit water uptake (Warren et al., 2015), therefore adding root biomass would merely add a greater respiratory cost to the plants and lead to lower allocation to photosynthetic material (leaves) and reduced net primary productivity (Sitch et al., 2003).

Over the relatively short (8-year) time series, there was no clear trend in AGC but the year-to-year anomalies (IAV) were mostly consistent in terms of sign between models and L-VOD AGC. The positive anomaly in 2011 was largely driven by the influence of the La Niña in Australia (Poulter et al., 2014; Haverd et al., 2016). These results could imply that LSMs capture the response to key drivers on carbon inputs via photosynthesis but not necessarily the allocation between above and belowground biomass or the processes that release carbon from these environments (e.g., respiration, fire, photodegradation). The L-VOD AGC anomalies have considerably higher amplitude than those of the LSMs. It is possible that the amplitudes of L-VOD-inferred anomalies may be exaggerated by variations in plant water stress, even though the L-VOD was averaged annually. The difference in amplitudes is furthermore sensitive to the biomass-transfer function used to calibrate the L-VOD data [e.g., linear, see this study (Supplementary Figure S3) and Brandt et al. (2018), or sinusoid, Fan et al. (2019)] and L-VOD data filtering. As L-VOD is a relatively recent product, methods for processing and analysing the data to best understand ecological functions are still being developed. Although longer time-series of VOD products exist, based on C, X, and Ku band microwave data, they have limitations such as increased sensitivity to water stress of foliage (Momen et al., 2017) and the necessity for intercalibration of different satellite instruments and measurements (Moesinger et al., 2020).

Simulated dryland vegetation carbon stocks varied between LSMs and decreased on average between 1901 and 2018 (Figures 6, 7) because reductions mainly due to land-use change exceeded the gains caused by CO₂ fertilization, while climate change had minimal effect in most models (Supplementary Figure S13). Dryland soil carbon stocks increased since 1970 despite being negatively impacted by land-use change in drylands such as the conversion of native woodlands to pasture and cropland which are often associated with degradation and can lead to decreased carbon input into soils (Cowie et al., 2011) (Supplementary Figure S14). These reductions in soil carbon from land-use change were counteracted by inputs from increased vegetation productivity due to CO₂ fertilization (Figures 6, 7, and Supplementary Figure S15) (Walker et al., 2021), resulting in soil carbon dominating the increasing trend in LSM simulated dryland ecosystem carbon storage (Figure 6). However, current ecosystem models are considered to poorly capture the complex interactions between CO₂ fertilization and soil organic carbon stocks (Terrer et al., 2021). With a lack of long-term observations, it is difficult to be confident in the role of dryland soils in climate-carbon feedbacks. This highlights the need for longer-term observational efforts to be able to verify these model estimates. Furthermore, we found large differences in the change in modelled carbon stocks of dryland ecosystems over centennial timescales between the different LSMs (Figures 6, 7). This level of disagreement in the change in carbon stocks when hindcast using relatively well-constrained climatology, CO₂ concentrations and land use is concerning because it undermines confidence in prognostic applications of these LSMs to even more uncertain future scenarios.

Fire is a critical process in many, but not all, dryland ecosystems, responsible for maintaining a stable ecological state in savanna ecosystems and is a major cause of carbon release (Bond et al., 2005; Lasslop et al., 2020). Land surface models increasingly include explicit representation of fire, including six of the twelve models considered here. Our analysis of these simulations found that the models with explicit fire did not show systematically better agreement with either GPP or AGC retrieved from remote sensing (Table 1). We found little systematic relationship between observed burn frequency and GPP residuals, although some models (JULES and CLASS-CTEM, ORCHIDEE-CNP, and OCN) did overestimate GPP in more frequently burned areas (Supplementary Figure S11). There was generally no relation between observed burn frequency and AGC residuals, apart from JULES and CLASS-CTEM which overestimated biomass in more frequently burned regions (Supplementary Figure S12). These results are consistent with the fire model intercomparisons (FireMIP), which found that while explicit fire substantially improved model-data correspondence in some regions (especially in functional drylands in South America that are excluded by our climatic definition of drylands) it worsened correspondence in other regions (Hantson et al., 2020). Vegetation models which lack explicit representation of fire contain compensatory biases which partially account for the effects of fire (Burton et al., 2019; Rabin et al., 2017). “LSMs that represent fires still do so insufficiently. For example, fire-enabled LSMs were not able to capture the global trend in burnt area (Andela et al., 2017). Uncertainties also remain regarding remote sensing observations with commonly used 500 m spatial resolution fire products significantly underestimating area burnt and fire carbon emissions (Ramo et al., 2021). Improvements to the representation of fire in LSMs are both needed and anticipated as better remotely-sensed fire products become available and we learn more about the compensatory biases present within current LSMs (Burton et al., 2019; Hantson et al., 2020; Lasslop et al., 2020).”

Simulations of terrestrial ecosystem processes in these LSMs struggle to capture many aspects of carbon dynamics in ecosystems that are subjected to severe water stress (Harper et al., 2020; MacBean et al., 2021). Limitations in model simulations of plant responses to water availability have been highlighted by MacBean et al. (2021) and Harper et al. (2020). For example, the plant functional types simulated in most models do not represent the ecophysiological adaptations present in most dryland vegetation communities (such as stomatal control and drought phenology), and simulated soils and rooting schemes are usually too shallow for these ecosystems. Furthermore, most LSMs have poor or no representation of biological soil crusts which are increasingly recognised to play critical roles in drylands including contributing ∼0.6 Pg yr⁻¹ of C to global net primary productivity (Chamizo et al., 2012; Rodriguez-Caballero et al., 2018). This is likely a reflection of a historical emphasis/bias of models to temperate forest phenology combined with a lack of empirical data. Differences in simulated vegetation dynamics also lead to erroneous estimates of fractional cover in key areas.

There are other processes thought to be important in drylands that LSMs do not represent. For instance, solar radiation interacting with plant litter causes photodegradation that emits gaseous carbon. Empirical evidence from decomposition (Austin & Vivanco, 2006; Day et al., 2018, 2019; Méndez et al., 2019; Berenstecher et al., 2020) and eddy covariance (Rutledge et al., 2010; Adair et al., 2017) studies suggests that photodegradation accounts for somewhere between 10% and 50% of gaseous carbon emissions in the semiarid ecosystems that dominate IAV and trends in dryland carbon dynamics (Poulter et al., 2014; Ahlström et al., 2015). Photodegradation is therefore an important aspect of carbon emission in seasonally arid ecosystems that is not currently represented in LSMs, despite their inclusion of insolation and litter pools. As the simulated long-term trends in ecosystem carbon stocks are dominated by cSoil (including litter) (Figures 6, 7), explicit representation of photodegradation would improve simulations of carbon dynamics in drylands.