Northern Hemisphere Snow Drought in Earth System Model Simulations and ERA5‐Land Data in 1980–2014

Low snow levels over the past few decades and predictions of a low‐to‐no snow future have spurred research into snow droughts, which pose a threat to water security and management. Systematic data‐model comparisons of snow drought have been lacking, hindering our understanding of the drivers of snow drought in the past. To address this gap, we analyzed snow drought events using standardized snow water equivalent index derived from monthly results of four numerical experiments using the E3SM Land Model (ELM) and ERA5‐Land data during the period of 1980–2014. Additionally, we compared snow drought duration calculated from models with those from the ERA5‐Land data during selected El Niño‐Southern Oscillation (ENSO) years. The numerical experiments were conducted with ELM driven by two prescribed atmospheric forcings, and with the coupled land‐atmosphere configuration of E3SM with and without plant hydraulics scheme feedback. Analysis reveals that 20%–30% of snow droughts occur due to factors other than above‐normal temperature and low snowfall, such as low soil moisture, warm soil temperature, and low relative humidity, etc., especially in high latitudes (50° North). Furthermore, our study highlights the exacerbating effect of ENSO events on snow drought conditions in various regions, despite some discrepancies between model and ERA5‐Land results. We also identified limitations of the coupled land‐atmosphere models in our current configuration in capturing the spatial patterns of snow droughts. This study underscores the challenge of predicting and mitigating snow drought and the need for a comprehensive understanding of the factors contributing to snow drought.


10.1029/2023JD039308
2 of 19 freshwater for many regions.Snow storage can affect ground temperatures, which is important for the analysis of geothermal resources assessment and use (Xu & Spitler, 2014), as well as for developing hydropower operation strategy (Mo et al., 2018).In addition to water supply, snow plays a critical role in global precipitation patterns (Field & Heymsfield, 2015).It also influences species distribution in ecosystems characterized by seasonal snow cover (Niittynen & Luoto, 2018).
Lack of snow accumulation, or snow drought, can have detrimental impacts, creating challenges for water management (Harpold et al., 2017).In addition, snow drought can present risk for food security and harm to wildlife and ecosystems (Barsugli et al., 2020;Keyser et al., 2023).For instance, the exceptional snow drought in the Italian Alps in the winter of 2021/2022 significantly reduced the energy potential of the hydroelectric reservoirs in Northern Italy and caused the deficit of melt water runoff (Koehler et al., 2022).Furthermore, snow drought can have a significant impact on forest phenology.Lou et al. (2021) found that snow drought advanced the start of growing season of forest by 6 days in a mid-high latitude area.In addition, snow drought can have severe implications for agriculture.For example, the record high snow drought in 2017/2018 in Afghanistan caused crop failures and devastated its agricultural sector (Huning & Aghakouchak, 2020).
According to the Glossary of Meteorology of American Meteorological Society, snow drought can be classified as "dry snow drought" due to the below-normal cold-season precipitation or "warm snow drought" due to a lack of snow accumulation despite near-normal precipitation.Warm snow droughts usually occur when warm temperatures prevent precipitation from falling as snow or result in unusually early snowmelt (Harpold & Kohler, 2017).However, it is important to note that these definitions may miss important components of how snow droughts originate and develop (Hatchett & McEvoy, 2018).The co-occurrence of warm and dry conditions is expected to increase in the future due to climate change.Studies indicate that dry snow drought is more likely to be followed by heatwave due to intensified soil drought and atmospheric aridity (X.F. Li & Wang, 2022).Additionally, snow droughts caused by both warm and dry winter conditions can lead to the most severe summer streamflow drought conditions (Dierauer et al., 2018).
Understanding the timing, location, and severity of snow drought is crucial for effective water resource management, mitigating the impacts of climate change, and achieving sustainable development goals related to water (Crausbay et al., 2020;Dierauer et al., 2021;Udall & Overpeck, 2017).Large-scale analysis and characterizations of snow droughts have been scarce until recently (Colombo et al., 2022;Cowherd et al., 2023;Hatchett et al., 2022;Huning & Aghakouchak, 2020).However, there has been a lack of systematic analysis of simulated results using Earth system models, which are essential for understanding the impacts of climate change on water availability and drought (Heavens et al., 2013).In this study, we aim to fill the gap by performing snow drought analysis on model results from an Earth system model.We will evaluate how the plant response to water availability affects snow drought in Earth system models as the water stress scheme (soil moisture scheme versus plant hydraulics scheme [PHS]) in the model can influence snow drought patterns through its impact on evapotranspiration and surface temperature.Compared to the soil moisture scheme, PHSs can provide a mechanistic linkage between plant transpiration and atmospheric water demand and soil water availability.
The amount of water released when a snowpack melts instantaneously, known as snow water equivalent (SWE), is commonly used to monitor snow drought.Various indices exist in the literature to classify snow droughts, such as the method proposed by Heldmyer et al. (2023) that clusters snow drought years into "warm," "dry," and "warm dry" categories using peak SWE and the ratio of basin-wide modeled peak SWE to accumulated (onset to peak) precipitation (SWE/P).Another method is using the US Drought Monitor "D" scale (Svoboda et al., 2002) and the standardized SWE index (SSWEI) (Huning & Aghakouchak, 2020).In this study, we adopt the latter method, which will be briefly introduced in the method section.
We compared our model results against a reanalysis data product to evaluate model skills and identify regions where the model needs improvement to accurately simulate snow droughts.We are particularly interested in snow droughts in El Niño-Southern Oscillation (ENSO) years, due to the lack of consensus on whether ENSO induces or relieves drought in snowpack conditions.For example, Thakur et al. (2020) identified that negative trends in western US from 1961 to 2018 SWE are largely associated with the ENSO phases, while another study by Tamaddun et al. (2017), covering a 56-year period , suggests an in-phase relationship between SWE and ENSO.

Snow Drought Characterization Method
We focused our study in the Northern Hemisphere.Following Huning and AghaKouchak (2020), we used 3-month SSWEI so the result is not sensitive to the occurrence of a single event.We considered grid where snow cover is greater than 0.5% and at least 75% of a given month for all years from 1980 to 2014 has a nonzero 3-month integration of the SWE value, which is given by where SWE m is the SWE value in month m from the model results.A m is the 3-month integrated SWE.m − 1 and m − 2 are the two previous months of the successive month m.A m s are the computed for November-April.
The probabilities associated with the 3-month SWE values in a given month m for all years were determined by computing the empirical Gringorten plotting position (Gringorten, 1963) of ordered SWE values (from smallest to largest): where m is the month, i is the rank of the nonzero SWE value, and N is the sample size.
SSWEI is then computed by transforming the empirical probability to the standard normal distribution as follows: where ϕ −1 is the inverse standard normal distribution.

Energy Exascale Earth System Model
Energy Exascale Earth System Model (E3SM) couples model components for atmosphere, land, ocean, sea ice, and river (Golaz et al., 2019;Leung et al., 2020).The atmospheric component of E3SM or EAM is developed based on the Community Atmosphere Model version 5 and uses the spectral element dycore with 72 vertical levels (Rasch et al., 2019).E3SM Land Model (ELM) is the land component of E3SM that branched from Community Land Model CLM4.5.It has a one-dimensional multilayer snow model.The number of vertical snow layers can be up to five depending on the snowpack depth.The snow processes in ELM include snow accumulation, melting, aging, densification, metamorphism, percolation and refreezing of water, and canopy snow interception (Hao et al., 2023;Toure et al., 2016).Each snow layer is characterized by its thickness, ice mass, liquid water content, temperature, and effective grain radius, and water can transfer between layers.A snow layer is initiated when the snow depth meets the minimum requirement of 0.01 m, and is subject to snow accumulation and melting.Excess water in a soil layer is allowed to percolate to the underlying layer when liquid water from rain or snowmelt reached its water-holding capacity (Toure et al., 2016).Snow cover fraction in ELM is estimated using the SWE and standard deviation of elevation (Hao et al., 2023), accounting for the parameterization that reflects the hysteresis in fractional snow cover for a given snow depth between accumulation and melt phases (Swenson & Lawrence, 2012).
In ELM, there are two approaches to regulate stomatal conductance when plants are under soil water stress (SMS).One is based on the soil moisture stress (SMS) formulation and the other is based on leaf water potential (ψ leaf ) formulation where ψ leaf is calculated from a PHS using Darcy's law to approximate the steady-state flow of water through the soil-plant-atmosphere continuum (Kennedy et al., 2019).
In SMS, the regulation of stomatal conductance is through the following stress factor β w,SMS : where ψ S,i is the soil water matric potential of soil layer i (m), ψ C and ψ O are the soil water potential (m) when stomata are fully closed and fully opened, respectively, and r i is the root fraction in soil layer i.When vegetation is not water stressed β is 1.0.
The stress factor (β w,PHS ) based on ψ leaf is modeled with the following sigmoidal function: ) (5) where P 50 is the water potential at 50% loss of stomatal conductance and c k is a shape-fitting parameter.
Note that β w,SMS and β w,PHS are also applied to reduce the maximum rate of carboxylation (V cmax ) in the model.
In coupled land-atmosphere simulations, differences in evapotranspiration due to different plant water stress representations may influence turbulence mixing in the planetary boundary layer, cloud formation, shortwave and longwave radiation, and atmospheric circulation.These factors may collectively affect the surface temperature and precipitation patterns, consequently leading to large differences in SWE and snow drought patterns.

Simulation Experiments
We performed four numerical experiments to simulate SWE, with the aim of evaluating the impact of climate forcing uncertainty on snow drought predictions.Two experiments were offline land simulations using ELM and global meteorological forcing data from observations, and two were global coupled land-atmosphere simulations using E3SM.The motivations for using forcings from both offline and coupled models are twofold: (a) to identify areas of agreement and disagreement between the model experiments and improve our understanding of the strengths and weaknesses of each model configuration; (b) to better assess the impacts of climate change on snow drought by taking into account the complex feedbacks and interactions between the land surface and the atmosphere.The first offline experiment utilized prescribed meteorological forcing data from phase 3 of the Global Soil Wetness Project (GSWP3) (Dirmeyer et al., 2006).GSWP3 has been widely used as a meteorological forcing data set in several climate impact assessments and snow modeling studies (Menard et al., 2019;Mengel et al., 2021;Muller Schmied et al., 2016;Schewe et al., 2019).The second offline experiment utilized GSWP3-W5E5, a combination of GSWP3 and W5E5, version 2 of Watch Forcing Data for ERA5 (WFDE5) (Cucchi et al., 2020;Lange et al., 2021).W5E5 precipitation is derived from WFDE5 rainfall and snowfall, bias-adjusted using monthly precipitation of version 2.3 of the Global Precipitation Climatology Project (Lange, 2019).GSWP3-W5E5 is considered to have a better representation of reality than GSWP3.Both forcing data have a 0.5° horizontal resolution, and the temporal resolution for GSWP3 and GSWP3-W5E5 are 3-hourly and 6-hourly, respectively.To conduct these experiments, we used ELM with a resolution of 0.5°, which was spun up for 80 years starting from 1900.We refer to these experiments as GSWP and W5E5, respectively.
We also conducted two global coupled land-atmosphere simulations using E3SM to investigate the impact of different plant water stress representations on model predictions.One of the simulations utilized an empirical SMS formulation, while the other used a formulation based on leaf water potential calculated from PHS.The atmosphere model was applied at a resolution of 1°, while the land model was applied at a resolution of 0.5°.In both simulations, sea surface temperature and sea ice were prescribed based on the observed climatology, while greenhouse gases, aerosols, and land use land cover were prescribed based on the transient historical forcing from the Coupled Model Intercomparison Project Phase 6 (CMIP6) input4mips data (Durack et al., 2018).For these simulations, the land model (ELM) was initialized using a restart file for 1980 from a previous simulation using the same transient model configuration.The atmosphere model was initialized using the atmospheric state of 1980 from an AMIP simulation driven by the same CMIP6 input4mips forcing data (Golaz et al., 2019).We refer to these experiments with an empirical water stress formulation and plant hydraulics as SMS and PHS, respectively.
The simulation results from all four experiments, spanning the period from 1980 to 2014, were evaluated against the gridded SSWEI calculated from ERA5-Land to evaluate the accuracy of meteorological forcing data for use in snow modeling and climate impact assessments and to identify areas where the E3SM model can be improved by testing different configurations to more accurately simulate hydrological processes such as snow droughts.

Benchmark Gridded SSWEI From ERA5-Land
To evaluate our model results, we use the latest land component of the reanalysis product ERA5-Land produced by the European Center for Medium-Range Weather Forecasts as a reference.This product is an enhanced global data set for ERA5 (Muñoz-Sabater, 2019; Muñoz-Sabater et al., 2021), which has been found to perform well 10.1029/2023JD039308 5 of 19 in the evaluation of nine gridded Northern Hemisphere snow SWE products as part of the European Space Agency Satellite Snow Product Intercomparison and Evaluation Exercise (SnowPEx) across a range of snow conditions (Mortimer et al., 2020;Shao et al., 2022).Although Mortimer et al. (2020) did not explicitly endorse any single product, they did note that ERA5 is one of the products that perform best across the range of snow conditions captured by the validation data set.ERA5-Land generally outperforms ERA5 over the US, while ERA5 performs better over Europe.It is noteworthy that the ERA5-Land SWE has good performance for heights ranging between ∼1,500 and ∼3,000 m (Muñoz-Sabater et al., 2021) and the best performance in the elevation interval of 300-400 m (Shao et al., 2022) but it is also known to have too large SWE sum estimates over the mountainous regions.
ERA5-Land is a reanalysis data set that provides the evolution of global land variables describing the water and energy cycles from 1950 to the present (Muñoz-Sabater et al., 2021).The global data set has a spatial resolution of 9 km.We downloaded monthly averaged data set of total precipitation, snowfall, snow cover, 2 m temperature in formats of General Regularly-distributed Information in Binary form from 1980 to 2014 via the user interface.
To ensure consistency in the spatial resolution of the modeled data, we interpolated the ERA5-Land data to a spatial resolution of 0.5°.

Data Analysis
During the period from 1980 to 2014, ENSO years were chosen for analysis, with 5 El Niño years having Oceanic Niño Index (ONI) greater than 0.7 and 5 La Niña years with ONI less than −1.5.The analysis focused on the Northern hemisphere, and the model-generated drought events and duration based on SSWEI were compared against that obtained from ERA5-Land.The state of temperature (cold "C," near normal "N," warm "W") and snowfall (dry "D," normal "N," wet "W"), or environmental conditions, under which each drought event occurs were analyzed by calculating the standardized index of temperature (STI) and snowfall (SSNI) similar to that of SSWEI.There are a total of nine combinations of temperature and snowfall conditions, namely, "CD" for cold and dry months when STI < −0.5 and SSNI < −0.5, "WD" for warm and dry months when STI > 0.5 and SSNI < −0.5, "WW" for warm and wet months when STI > 0.5 and SSNI > 0.5, "NN" for near normal months when STI and SSNI are both between −0.5 and 0.5, and so on.These combinations of temperature and snowfall conditions were considered to identify droughts caused by warm temperature, low snowfall, or both warm and low snowfall.If a snow drought event falls into the category that is neither warm nor dry, the standardized indices of soil moisture and soil temperature of soil layer 1 (0-7 cm) from ERA5-Land and the top 10 cm soil layer from E3SM will be calculated to further evaluate the driver of snow droughts.The durations of snow droughts in months were calculated for each product for the selected ENSO years to identify regions that are prone to snow droughts lasting more than two consecutive months or hotspots of snow droughts.Statistical significance (P < 0.05) between models and ERA5 for environmental conditions that contribute to snow droughts and drought durations were conducted with one-way analysis of variance (ANOVA) using Python function scipy.stats.f_oneway.Post hoc comparisons between groups were conducted with Tukey's honestly significant difference (HSD) test using python function statsmodels.stats.multicomp.tukeyhsd.

Model Bias
Our analysis of long-term climatology reveals a negative bias between the results obtained from offline land simulations using prescribed atmospheric forcing (i.e., GSWP and W5E5) and ERA5 in mountainous areas compared to the rest of the land (Figure 1).This is mostly due to warmer temperature and lower precipitation in GSWP and W5E5 in those regions compared to those in ERA5, as shown in Figures S1 and S2 in Supporting Information S1.
The bias in the western US is improved by using W5E5 as forcing.The large SWE biases between ERA5 and the ELM simulation driven by GSWP3-W5E5 in eastern US and China are due to the large difference in long-term average precipitation between W5E5 and ERA5 in those regions, as depicted in Figure 6c in Cucchi et al. (2020).Except for mountainous areas, the SWE bias from the coupled models is generally low in high latitudes, but high in low latitudes (as shown in Figure 1).Note that the legend in Figure 1 is on a log10 scale.The actual magnitude of SWE differences in the low latitudes is generally lower compared to the high latitudes.The bias is consistent between the two coupled model configurations.It is noteworthy that SWE bias from the coupled models is generally smaller than from the offline land models in latitudes 40°N and above. 10.1029/2023JD039308 6 of 19 The drought classification based on SSWEI from the offline land models and ERA5 exhibit certain spatial similarities regarding the occurrence of snow droughts and the progression of droughts over time.However, there are variations in the severity of droughts between these models.For example, in Figure 2 depicting the year 1998 (a strong El Niño year), it is evident that the droughts from the offline land models in the US were more severe compared to ERA5.Among the different models, GSWP demonstrates the closest agreement with ERA5 in terms of drought severity.
On the other hand, the spatial distribution of droughts from the coupled models differs significantly not only from each other but also from ERA5, as can be seen in the Eurasian region in 1998 (Figure 2).These differences highlight the limitations of the coupled models in our current configuration in capturing the spatial patterns of droughts, likely due to differences between the simulated atmospheric environment and the atmospheric forcing used in the offline simulations.

Environmental Conditions Contributing to Snow Drought
To gain a better understanding of the causes of each snow drought event, we compared the temperature and snowfall among the three reanalysis data sets (ERA5-Land, GSWP3, W5E5) and those from the coupled simulations using E3SM, initially assuming snowfall and temperature are the main drivers of snow droughts.The grids experiencing snow droughts from 1980 to 2014 were divided into nine categories based on the standardized temperature index (STI) and standardized snowfall index (SSI).The categories were labeled as "WD" (warm, dry) when STI > 0.5 and SSI < −0.5, "NN" (near normal) for STI and SSI between −0.5 and 0.5, "CW" (cold, wet) for STI < −0.5 and SSI > 0.5, and so on.Figure 3 displays the spatial conditions in March and April 2008 as an example.The analysis revealed that some instances of snow droughts occurred when the temperature was low and snowfall was above normal, that is, categories marked as cold, near normal snowfall ("CN"), near normal temperature, wet snowfall ("NW"), near normal temperature, near-normal snowfall ("NN"), and cold, wet ("CW") in Figure 3.For those instances that do not fall into the "dry" or "warm" or "dry and warm" category, we then evaluate the standardized indices of soil moisture and soil temperature for that pixel.To determine the fraction of snow drought under each of these environmental conditions for each pixel, we calculate how many snow droughts a specific condition occurred divided by the total number of snow droughts for that pixel.Note that many of these conditions can co-occur.We attribute the main driver of a snow drought following our analysis sequence.For instance, when both low snowfall and high soil temperatures coincide during a snow drought, our analysis attributes the snow drought to low snowfall.
Although the dominant conditions for snow drought are typically due to warm temperature and/or low snowfall, a significant fraction of snow droughts occurred under normal or below-normal temperature and normal or above-normal snowfall conditions, as shown in Figure 4. Droughts in the lower latitudes were mostly due to combined warm temperature and low snowfall or WD condition, while in the high latitudes, they were mainly caused by low snowfall or dry condition.While the models and ERA5 share some similarities in depicting snow droughts attributable to warm, dry, and WD conditions, close inspection reveals localized variations (Figure 4).Specifically, when examining dry, warm conditions, the spatial patterns show greater consistency among the three reanalysis data sets than with those from the coupled simulations.Under WD conditions, ERA5 exhibits a notably distinct spatial distribution of snow drought, particularly evident in the Tibetan Plateau region.These observations between models and ERA5 for these three types of snow drought are also confirmed by the probability density function (PDF) of snow drought under the warm, dry, and WD environments (Figure 5).The fraction of snow drought due to warm temperature is less than 0.2 (Figure 5a), whereas low snowfall results in a fraction greater than 0.3 (Figure 5b).When warm temperatures and low snowfall occur simultaneously, the fraction distribution exhibits a bimodal pattern, typically ranging between 0.15 at high latitudes and 0.5 at low latitudes (Figures 4 and 5c).On average, the fractions of snow drought attributed to warm temperatures are smaller, while the occurrences of reduced snowfall (dry) are more pronounced, comparing results from the reanalysis forcing to the results obtained from the coupled model simulations (see Table 1).
Further analysis of the standardized index of soil water and soil temperature for some fractions of snow droughts not explained by warm temperature or low snowfall revealed that soil moisture was below normal, the soil temperature was above normal, or both, for those snow droughts (Figures 4  and 5d-5f).The combined average fraction attributable to soil water and/ or temperature comprises approximately 20% in ERA5 and about 30% in model results (Table 1).These findings suggest that either warm or dry soil can increase the melting rate of the snowpack or the pre-existing snowpack is too low to recover from the drought.A portion of snow droughts was explained by below-normal relative humidity, which can increase the rate of sublimation and reduce the snowpack.However, there remains an unexplained fraction (∼7% on average; Table 1) for both ERA5 and the models (Figure 5g).The models show the best agreement with ERA5 in the fraction of snow droughts caused by warm soil or dry soil.For the models, the fraction of snow droughts that cannot be explained by temperature or snowfall in the month when snow drought occurred is typically around 0.3 and for ERA5, this is 0.16.The fraction distribution of snow droughts due to dry and warm soil conditions from ERA5 differs the most from that of the models, showing a narrower fraction range (Figure 5f), and the mean fraction is 0.07 compared to 0.14 from the models.This difference may be attributed to differences in how soil hydrology and surface runoff are represented in ERA5-Land and ELM, as these processes are not well constrained by observations or atmospheric forcing.
Since snowpack has a strong dependence on surface elevation, we categorized elevations into six distinct groups to evaluate how the drivers of snow drought may vary with elevation.The PDFs of the fraction of snow drought attributed to warm atmospheric conditions from the ERA5 data set exhibit a wider spread at elevations higher than 1,000 m, with the highest median occurring in the elevation bin of 3,000-6,000 m.Conversely, the median fraction attributed to dry atmospheric conditions decreases with elevations exceeding 600 m (see Figure 6a).The median fraction attributed to warm and dry atmospheric conditions is higher for high-elevation bins compared to low-elevation bins (see Figures 6b and 6c).Notably, the shape of the PDFs representing fractions attributed to other environmental conditions (except for dry soil) within the elevation bin of 3,000-6,000 m differs significantly from those of the other elevation bins (see Figures 6e-6g).
The fractions of snow drought attributed to different factors, excluding "others," in elevation bins below 1,000 m (Figure 6) exhibit similarities across various data set groups from model experiments.Specifically, GSWP demonstrates similarities with W5E5, while SMS aligns with PHS.In contrast to ERA5, the medians of fractions under dry atmospheric conditions from GSWP and W5E5 are the lowest in the highest elevation range (Figure 6b).Similar to ERA5, the medians are lower in elevation ranges above 1,000 m, albeit with a smaller magnitude (or slower rate of decrease).Like ERA5, the medians of fractions attributed to warm atmospheres in SMS and PHS are the highest in the elevation bin of 3,000-6,000 m, whereas GSWP and W5E5 exhibit smaller medians for this elevation group (Figure 6a).All data sets exhibit significantly higher median fractions attributed to warm and dry atmospheric conditions when the elevation exceeds 1,000 m (Figure 6c).

ENSO Snow Drought Signatures
In the remaining subsections, we will focus on snow drought identified in ENSO years to evaluate whether ENSO induces or relieves drought in snowpack conditions and in what regions.We calculated ENSO signatures for each model and ERA5 by subtracting the annual mean of precipitation, temperature, and SWE during selected ENSO years from the annual mean of the same variables during all simulation years.We then compared the model signatures of ENSO with the ERA5 signatures to assess the agreement between the models and the reference.In ERA5, El Niño is characterized by warmer temperature in Northern America, lower precipitation in the US, and reduced SWE primarily in the northern US (Figure 7).In the west of Central Asia, it corresponds to cooler temperatures and higher precipitation, while in eastern Russia reduced SWE corresponds to warmer temperature and higher precipitation.In La Niña years, a decrease in SWE primarily occurs in northern North America, which is associated with cooler temperatures and lower precipitation.In contrast, in northern Asia, a decrease in SWE is linked to mixed conditions of temperature and precipitation.El Niño has a greater impact on SWE reduction in North America than La Niña.
The magnitudes of signature differences in the W5E5 and GSWP models are similar, except for the more pronounced temperature signal observed in eastern Canada in the W5E5 model (Figure S3 in Supporting Information S1).The SMS model shows heavier precipitation in southwestern US, eastern North America, Europe, and eastern and southern China during El Niño is a more negative SWE bias compared to ERA5 (Figure S5 in Supporting Information S1), although some exceptions exist.For instance, in northern Europe, the negative SWE bias from SMS corresponds to a positive SSWEI bias, whereas negative SWE bias in northern US primarily corresponds to positive SSWEI bias in El Niño years, suggesting that the snowpack recovered quickly from snow drought in some events in the model.In other words, the low SWE events were sporadic and might not lead to snow droughts.SSWEI exhibits a larger negative bias from SMS compared to the other models, particularly in La Niña years.In Canada, PHS shows a larger positive SSWEI bias in El Niño years.

ENSO Effect on Snow Drought Duration
To investigate the lasting impact of ENSO on snowpack, we identified consecutive months classified as snow droughts (SSWEI < −0.5) and calculated the duration of snow drought for each year.The models and ERA5 agree that ENSO events generally prolong the snow droughts compared to long-term climatology (as shown in Figure 9), with some locations experiencing droughts that last more than 2 months longer.
The snow drought duration differences between ENSO years and long-term climatology reveal a wider spread of data points with a range of −2 to 4 (Figure S6 in Supporting Information S1).There is a concentration of data points around positive values of 0-2, indicating longer drought durations during ENSO events compared to the long-term climatology in some regions (Figure S6 in Supporting Information S1).The long, narrow shape in the positive range for ERA5 and the models indicates ENSO events worsen snow droughts.The results of the statistical analysis indicate that there are significant differences in snow drought duration between groups during ENSO events (1-way ANOVA, P < 0.05.Post hoc analysis using the Tukey HSD test was conducted to further explore the differences between groups. In El Niño years, the results showed that the mean drought duration from W3E5 was the longest (with a mean difference between ERA5 of 0.09), and PHS was the shortest (mean difference = −0.23)among the ERA5 and model results (Figure S7 in Supporting Information S1).There was no statistically significant difference between ERA5, GSWP, and SMS.In Europe, the mean drought duration for PHS was the longest (mean difference = 0.19), and SMS was the shortest (mean difference = −0.35).There was no statistically significant difference between GSWP and W5E5.In Asia, drought duration from GSWP was the longest (mean difference = 0.22), and SMS was the shortest (mean difference = −0.16).There was no statistically significant difference between GSWP and W5E5.
In La Niña years, the results showed that the mean drought duration for SMS in North America was the longest (mean difference = 0.4), and PHS was the   Each row corresponds to specific environmental conditions, including "warm atmosphere," "dry atmosphere," "warm and dry atmosphere," "dry soil," "warm soil," "dry and warm soil," and "others," listed from top to bottom.
shortest (mean difference = −0.15)among model and ERA5 results (Figure S8 in Supporting Information There was no statistically significant difference between ERA5 and GSWP.In Europe, the longest duration was observed in PHS (mean difference = 0.48) and the shortest was observed in GSWP (mean difference = −0.24).In Asia, the longest duration was observed in SMS (mean difference = 0.22), and the shortest was observed in GSWP (mean difference = −0.23).
During El Niño events, warmer and drier conditions typically occur in various regions, including parts of the US, Europe, and Asia.According to SSWEI from ERA5, snow droughts of longer duration can be found in the Cascade Range in the western US, parts of the eastern US, northwestern Canada, central and northern Europe, as well as northern Russia.Although drought durations from the models generally agree with ERA5-Land for most of the regions, they differ in terms of exact locations and durations.Conversely, La Niña can have the opposite effect in the regions mentioned above, leading to shorter durations of droughts, according to ERA5 (as shown in Figure 9).ERA5 also indicates more frequent droughts in Central and Eastern Asia during La Niña years.
The offline models generally agree better with ERA5 than the coupled models in terms of the spatial distribution of longer durations, except for the high mountain Asia region.The discrepancy of drought durations in the ENSO years between the GSWP3 and W5E5 experiments is small compared to the other experiments.W5E5 has a relatively larger area of droughts with longer durations in La Niña years in parts of Central Asia and matches ERA5 better.SMS simulates longer average durations in central North America during El Niño years, and in Canada and eastern Russia in La Niña years.PHS simulates shorter durations in Canada, and longer durations in northeastern Russia during El Niño years, as well as shorter durations in the Cascade Range in North America and northern Russia, and longer durations in Europe during La Niña years.

Environmental Conditions Contributing to Snow Droughts in ENSO Years
Based on the Tukey HSD test, during El Niño years, the mean fraction difference between ERA5 and W5E5 was not significant for all environmental conditions except for dry and warm soil.The mean fraction difference 10.1029/2023JD039308 13 of 19 between PHS and ERA5 was larger than that of the long-term climatology for almost all conditions except for dry and warm soil.Figure S7 in Supporting Information S1 compared to Figure 5 shows this difference.The mean fraction for each condition increased more than 0.1 compared to the long-term climatology, and the mean fraction of unknown contributors increased the most.
During La Niña years, the Tukey HSD test indicates that the mean fraction among ERA5, GSWP, and PHS was not significantly different under warm temperature, dry soil, warm soil, and unknown conditions.In comparison to the long-term climatology, data were spread out over a larger range of fractions during La Niña years, suggesting a more even distribution across the space (Figure S8 in Supporting Information S1 compared to Figure 5).The mean fraction for each condition also increased more than 0.1 compared to the long-term climatology, and mean fraction of unknown contributors increased the most.Compared to the El Niño years, the mean fractions contributed to low snowfall decreased from 0.44 to 0.38.
In El Niño years, the drought in North America is primarily driven by warm and dry atmospheric conditions, while low snowfall is the main driver in eastern Asia, and warm temperatures is the main driver in Europe, according to the results from ERA5 (as shown in the fraction difference between El Niño years and full simulation period in Figure S9 in Supporting Information S1).Warm and dry soil also contributes to the prolonged drought in Europe and eastern Asia (Figure 9), which is mostly supported by the models.The longer duration of drought  in North America from SMS compared to the other models is due to more occurrences of warm and dry atmospheric conditions.On the other hand, the shorter duration of drought in Canada from PHS is because of the low occurrences of warm or dry atmospheric conditions.
In La Niña years, according to ERA5, the long duration of drought in western Canada is mainly due to low snowfall and dry and warm soil, while warm temperatures and other causes contribute to drought in central Asia (Figure 9 and Figure S10 in Supporting Information S1).The snow drought in eastern China is primarily due to warm and dry atmospheric conditions according to ERA5.The models are generally in agreement with ERA5, with a few exceptions.SMS and PHS indicate that the long drought duration in western Canada was mainly due to frequent low snowfall.In eastern North America, the drought duration is longer due to the frequent occurrence of warm and dry atmospheric conditions.In northeastern China, the drought is mainly due to warm and dry atmospheric conditions based on results from ERA5, GSWP, and W5E5, while it is due to low snowfall from SMS, and warm and dry soil from PHS. PHS has a shorter drought duration in the Cascade regions in North America due to the low occurrences of warm or dry atmospheric conditions.

Discussions
Using E3SM, we performed four numerical experiments including two offline land simulations driven by two reanalysis atmospheric forcing data, and two global coupled land-atmosphere simulations with different plant water stress representations.The goals were to understand the drivers of snow drought, including contrasting the drivers during ENSO and long-term climatology, and to benchmark snow drought distribution from model results with those from the ERA5-Land data and explore the sources of uncertainty on snow drought simulations.

Factors Contributing to Snow Droughts
Snow drought is a complex phenomenon that can be caused by both direct and indirect factors.Direct factors include low snowfall, high temperatures, and wind events that can cause snow sublimation or drifting.Indirect factors including soil moisture, soil temperature, and atmospheric water demand also play important roles in snow drought.Generally, dry soil conditions and high soil temperature can result in increased snowmelt or reduced snow accumulation, both leading to reduced snowpack.Based on our analysis, we found 20%-30% of land fractions experiencing snow drought have sufficient snowfall and low temperatures (Figures 4 and 5), but the soil is dry or the soil temperature is high, making snow unable to accumulate and melt quickly, further exacerbating the occurrence of snow drought.Additionally, atmospheric water demand, driven by factors such as temperature, humidity, and wind, can indirectly affect snow drought by increasing the amount of water lost from the snowpack through sublimation.These direct and indirect factors can interact in complex ways, making it challenging to predict and mitigate snow drought.To fully understand and characterize snow drought for developing effective strategies for water management and climate adaptation, it is important to consider a range of factors that can influence the snowpack besides snowfall and temperature, and how well they are simulated by models.

Impact of ENSO Events on Snow Drought
Our findings indicate that the impact of ENSO events on snow drought varies depending on the location, with differences in the drought duration of aggravation and relief.In North America, both El Niño and La Niña events can induce prolonged snow drought, but the duration is less pronounced during La Niña events (Figure 9).In Russia, prolonged snow drought is more likely to occur during El Niño events, while central Asia is more prone to longer snow drought during La Niña events.The remaining land areas experience relief from snow drought during ENSO events.

Discrepancies in Snow Drought Between Models and ERA5
The majority of snow droughts are driven by atmospheric forcings.While there are some similarities in the spatial distribution of snow droughts between the models and ERA5, there are discrepancies in locations that are more prone to severe and prolonged snow droughts using different atmospheric forcings, especially in the high mountain areas.For example, when using the improved atmospheric forcing GSWP3-W5E5, ELM largely underpredicts snow drought in High Mountain Asia compared to ERA5.Noticeable differences between the 10.1029/2023JD039308 16 of 19 coupled simulations and offline simulations are expected due to differences in the atmospheric environments between the coupled simulations and the atmospheric forcing used to drive the offline simulations, as the former are free-running simulations constrained only by observed sea surface temperature and sea ice conditions.
It is worth noting that ERA5-Land's snowpack is determined solely by using the atmospheric forcings obtained from ERA5 to drive the H-TESSEL land surface model (Balsamo et al., 2009).Although no data assimilation was used in ERA5-Land, the atmospheric forcings were directly constrained by assimilated observations (Räisänen, 2021).However, ERA5-Land has been found to show especially large overestimation of SWE in the European Alps (Monteiro & Morin, 2023).On the other hand, Lei et al. (2022) found snow depth of ERA5-Land represents snow absence well and compares reasonably well with in situ observations of snow depth over the Tibetan Plateau.Uncertainties and biases in ERA5-Land suggest some uncertainties in interpreting the comparison of ELM with ERA5-Land and highlight the importance of validating model outputs with ground observations and considering the limitations and uncertainties of both models and reference data sets.The use of a single reference product, ERA5-Land, also presents a limitation in our evaluation process in terms of the model performance.In our future research endeavors, we will address this limitation by incorporating a broader array of SWE products to ensure a more comprehensive and robust assessment of our model results.
These discrepancies in model predictions may be caused by several factors, including inadequate model resolution to capture the complex terrain and topography-induced variability in meteorological conditions, simplified parameterization schemes that cannot accurately represent the physical snow processes, uncertainty in atmospheric forcing, and models not fully capturing the complex and nonlinear feedback interactions between snow pack, vegetation, and climate (Günther et al., 2019;Kim et al., 2021;Varga & Breuer, 2023).

Plant Hydraulic Scheme Can Change Snow Drought Patterns
In coupled land-atmosphere simulations, differences in plant response to climate forcing due to differences in how models simulate the plants response to water stress may affect the partitioning between the surface sensible and latent heat fluxes, which in turn can impact clouds and precipitation, modulating atmospheric latent heating and atmospheric circulation (Fang & Leung, 2022).These changes can further modify regional surface air temperature and humidity by influencing regional vapor pressure deficit.Green et al. (2017) provide evidence that feedbacks between vegetation and atmosphere can impact climatic conditions and explain up to 30% of variance in precipitation and surface radiation.Although there are similar long-term spatial patterns between model and ERA5 for total precipitation and temperature in both SMS and PHS (Figures S1 and S2 in Supporting Information S1), there is a notable difference in temperature, precipitation, and SWE during ENSO years (as shown in Figures S3-S5 in Supporting Information S1).This suggests that PHS can simulate the strong feedback between soil hydrology and root water uptake for transpiration, which in turn affects the regional air temperature, precipitation, and consequently dynamics of the snowpack.

Additional Factors Affecting Snow Drought
There are many factors we did not consider in this study that can also affect snow drought.Human activities perturbing land conditioning by deforestation, urbanization, and agricultural practices can alter the surface albedo and energy balance (He et al., 2014;Ouyang et al., 2022), which can affect the timing and amount of snowfall and snowmelt.Additionally, urbanization can increase the amount of impervious surfaces, changing the hydrological cycle, and resulting in reduced soil moisture and snowpack accumulation (Guo et al., 2021;Perryman & Dixon, 2013;Scheuer et al., 2017).Wildfires can also have a significant impact on snow drought (Kampf et al., 2022;Smoot & Gleason, 2021).Wildfires can burn off the vegetation cover and burned materials fall onto the snow surface, which can reduce the albedo, and increase the amount of shortwave radiation absorbed by the snowpack, leading to increased snowmelt rates, and decreased snowpack accumulation (Gleason et al., 2013).These human impact and wildfire highlight the importance of considering these factors in efforts to understand and mitigate snow drought.

Conclusions
Understanding the threat that snow droughts pose to water security and management is crucial, which is why the study of snow droughts is very important.In this research, we analyzed snow droughts using model results and compared them with reference data derived from ERA5-Land, which fills the gap of lacking global systematic data-model comparisons.Our study showed that ENSO events can aggravate snow droughts in various regions.Specifically, El Niño events were found to potentially prolong snow droughts in Russia, whereas La Niña events may contribute to extended snow droughts in Central Asia.The occurrences of snow drought are predominantly due to low snowfall, warm temperature, or both.However, about 20%-30% of snow droughts in our analysis can occur due to factors other than above-normal temperature and low snowfall, such as low soil moisture, warm soil temperature, low relative humidity.Therefore, it is important to take into account all the factors that contribute to drought conditions besides temperature and snowfall conditions.It is also important to consider the local conditions and how they interact with larger climate patterns associated with modes of climate variability such as ENSO to fully understand the causes of drought.Furthermore, different factors may respond differently to climate change, and modulate snow drought differently in the future.Future investigations considering all the factors contributing to snow drought can inform better management strategies for water resources and ecosystem services in the face of increasing snow drought risk.
Although snow droughts are generally well simulated by the models in terms of the broad spatial patterns, this study also revealed a significant discrepancy between the models and ERA5-Land.While E3SM accounts for land-atmosphere interactions, there is a need for improvement in how it simulates snow processes, assuming ERA5-Land SWE is accurate.Offline land simulations, even with improved atmospheric forcing, still struggle to predict snow droughts in High Mountain Asia because there is not enough data that meet the criteria to calculate SSWEI.The simulated SWE from offline models in this region is underrepresented compared to ERA5-Land.These results suggest that the model may poorly represent snowpack dynamics in the area, which needs to be further explored with more observations in the future.
In conclusion, improvements in atmospheric forcing and land surface models can help in accurately simulating and predicting snow droughts, which in turn can aid in decision-making for water resource management.Further research is needed to develop better approaches for modeling snow processes and soil moisture dynamics in Earth system models.

Figure 3 .
Figure 3.The environmental conditions during snow droughts in March (left panels) and April (right panels) 2008.Each row represents analysis from different products.In the caption, the left letter represents temperature condition (cold, normal, and warm) and the right represents snowfall condition (dry, normal, and wet).

Figure 5 .
Figure 5. Probability density function of snow drought fraction under different environmental conditions from 1980 to 2014.

Figure 6 .
Figure 6.Probability density function of the snow drought fraction within different elevation bins and under various environmental conditions from 1980 to 2014.The width of each violin corresponds to the data density.The horizontal dashed lines from the bottom in each violin represent the first quartile, median, and third quartile, respectively.Each row corresponds to specific environmental conditions, including "warm atmosphere," "dry atmosphere," "warm and dry atmosphere," "dry soil," "warm soil," "dry and warm soil," and "others," listed from top to bottom.

Figure 7 .
Figure 7. El Niño-Southern Oscillation (ENSO) signatures from ERA5 (annual mean in ENSO years-annual mean in all simulation years).

Figure 8 .
Figure 8. Difference in El Niño-Southern Oscillation signature of standardized snow water equivalent index between ERA5 and models.

Figure 9 .
Figure 9. Spatial distribution of the difference in average snow drought duration (months) over selected El Niño-Southern Oscillation years and long-term climatology (El Niño on the left, and La Niña on the right).