Document Type : Research Article
Authors
Department of Physical Geography, Faculty of Geographical Sciences, Kharazmi University, Tehran, Iran
Abstract
Snow water equivalent (SWE), snow depth (SD) and bulk snow density (ρs) are the fundamental state variables of snow hydrology, yet a simultaneous, country-wide assessment of all three has rarely been undertaken for Iran. This study analyses their spatio-temporal variability using monthly ERA5-Land reanalysis fields at ~0.1° (~9 km) resolution for 1980–2024 (15,547 land grid cells). Bivariate relationships were evaluated with scatter plots, regression and Pearson and Spearman coefficients computed for four sample definitions, taking the definitional dependence SWE = SD × ρs/1000 into account; principal component analysis (PCA) and K-means clustering were then applied, with the number of clusters assessed by the elbow, silhouette and Calinski–Harabasz criteria and a two-variable sensitivity test. Temporal change was examined with the modified Mann–Kendall test, Sen's slope and the Pettitt test on spatially averaged series, and with cell-wise trend maps and two-period (water years 1981–2003 vs 2004–2024) difference maps. Snow storage is extremely concentrated: 27.3% of Iran has a long-term mean SWE above 1 mm, half of the storage lies in 3.7% of the area (Gini = 0.87), and the maximum (0.155 m w.e.) occurs in the Zardkuh–Dena sector of the central Zagros. Depth and SWE are almost collinear (r = 0.99 for long-term cell means), partly by construction, whereas density is positively but nonlinearly related to depth and its correlations depend strongly on the sample (r = 0.62–0.87). PC1 explains 90.8% of the variance, and four snow classes were distinguished along a continuous gradient; the 184 cells of the mountain-snowpack class (1.2% of the area) hold 25.8% of the storage. December–March mean SWE declined significantly (Z = −3.73; Sen's slope −1.9 mm w.e. per decade), with a change point after 2008 that halved the winter mean; April snow declined significantly, whereas November showed no trend. Significant winter SWE decreases cover 80.3% of the snow-affected area and no cell shows a significant increase. Because the results derive from a single reanalysis product, they describe ERA5-Land and require validation against in-situ observations.
Extended Abstract
Introduction
Snow water equivalent (SWE), snow depth (SD) and bulk snow density (ρs) are the three fundamental state variables of snow hydrology. SWE, the product of depth and density, determines the volume of water stored in the snowpack and released during the melt season, whereas density is the key parameter for converting the far more abundant depth observations into SWE (Sturm et al., 2010; Hill et al., 2019; Che, Dai & Li, 2025). Iran is a predominantly semi-arid country whose agricultural, domestic and industrial water supply depends heavily on the seasonal snowpack of the Alborz and Zagros ranges and the northwestern highlands, and mounting evidence points to declining snow depth, snow-cover duration and snowmelt-derived water resources across the country (Nouri & Homaee, 2021; Sadeqi et al., 2024). Iranian snow research, however, remains dominated by basin-scale or station-based case studies; the only country-wide assessment of all three variables (Vafakhah, 2018) relied on nine stations, and independent studies of snow density are scarce. The international literature is also inconsistent regarding the depth–density relationship, which has been reported as negative at short time scales (McCreight & Small, 2014), absent in the deep Alpine snowpack (Guyennon, Valt, Salerno, Petrangeli & Romano, 2019) and positive for shallow snow in a dataset of more than 6,000 Northern Hemisphere stations (Che et al., 2025). In addition, ERA5-Land has reproduced winter snow-depth variations over northwestern Iran better than MERRA-2 (Majidi Karhroudi, Gharaylou & Sabetghadam, 2024; Salahi & Shahi, 2026); nevertheless, it has not yet been used to examine all three snow variables jointly at the national scale.
The present study therefore pursues three objectives: (i) to describe the spatial structure and the degree of spatial concentration of SWE, SD and ρs across the whole of Iran; (ii) to quantify the interrelationships of the three variables while explicitly accounting for their definitional dependence (SWE = SD × ρs/1000) and for the effect of the sampling unit on the correlation coefficients; and (iii) to test for monotonic trends, change points and spatial patterns of change during 1980–2024. The novelty of the study lies in four aspects. First, it is the first country-wide, grid-based (~9 km) assessment of all three snow variables in Iran, covering 15,547 land cells rather than a few stations or a single basin. Second, it treats the definitional dependence among the variables explicitly, so that structurally inflated correlations are not presented as physical findings. Third, by recomputing the coefficients for four sampling units, it shows that part of the apparent inconsistency in the published depth–density relationship can be explained by differences in scale and sampling unit. Fourth, it combines measures of spatial concentration (the Lorenz curve and the Gini coefficient) with multivariate classification and with trend and change-point analysis in order to identify where Iran's snow storage is concentrated and where it has declined most.
Material and Methods
The study domain is Iran (≈1,648,195 km²; 25°–40°N, 44°–63.5°E). Monthly-mean ERA5-Land fields (Muñoz-Sabater et al., 2021) of snow depth (sde, m), snow depth water equivalent (sd, m w.e.) and snow density (rsn, kg m⁻³) were obtained from the Copernicus Climate Data Store at a resolution of 0.1° (~9 km) for January 1980–December 2024 (540 months); the incomplete year 2025 was excluded. Density was taken directly from rsn, and the relation SD ≈ SWE × 1000/ρs was satisfied with a median ratio of 0.997. Quality control covered 15,547 land cells (8,395,380 values per variable). No missing or physically invalid values were found, very small negative values resulting from the numerical storage of zeros were set to zero, and 68.1% of all density values equal the 100 kg m⁻³ minimum of the single-layer snow scheme (Dutra et al., 2010), which is assigned to snow-free time steps.
The analysis proceeded in five stages. (1) Bivariate relationships were assessed on long-term cell means with scatter plots, a saturating exponential density–depth fit and Pearson (r) and Spearman (rs) coefficients. Because SWE is the product of depth and density, part of the correlation among the variables is structural (Kronmal, 1993); the coefficients were therefore recomputed for four sampling units (long-term cell means, monthly spatial means, all cell-months and snow-present cell-months) and are reported descriptively, since neighbouring cells are spatially autocorrelated. (2) Principal component analysis (PCA) was applied to the standardized variables, and K-means clustering was performed for K = 2–8, with the number of clusters evaluated by the elbow, silhouette (Rousseeuw, 1987) and Calinski–Harabasz (Caliński & Harabasz, 1974) criteria; a sensitivity analysis repeated the procedure with SD and ρs only, and the agreement between the two partitions was measured with the adjusted Rand index (ARI; Hubert & Arabie, 1985). (3) The spatial concentration of storage was quantified with a Lorenz curve and the Gini coefficient. (4) Spatially averaged December–March, April, November and annual series of SD and SWE were tested with the Mann–Kendall test and its modification for autocorrelation (Hamed & Rao, 1998), Sen's slope (Sen, 1968) and the Pettitt change-point test (Pettitt, 1979) at α = 0.05. (5) Cell-wise Mann–Kendall and Sen's slope maps of December–March SWE and of each calendar month were produced for the 5,695 snow-affected cells (mean winter SWE ≥ 1 mm), together with difference maps between water years 1981–2003 and 2004–2024, the division used in the SWEI snow-drought analysis of Iran (Hajivand Paydari, Saligheh, Alijani & Karbalaee, 2025).
Results and Discussion
Snow storage is strongly orographic and extremely concentrated. Only 27.3% of Iran has a long-term mean SWE above 1 mm; half of the national storage lies within 3.7% of the area, 80% within 10.6% and 90% within 16.2% (Gini = 0.87). The absolute maximum (0.155 m w.e.) occurs at 32.5°N, 50.0°E in the Zardkuh–Dena sector of the central Zagros, the headwaters of the Karun and Zayandeh-Rud rivers, with an elongated secondary centre along the central Alborz and a third centre in the northwestern highlands (Sahand–Sabalan).
Depth and SWE are almost collinear (r = 0.99, rs = 1.00 for cell means), which is largely structural. Density increases nonlinearly with depth (r = 0.81, rs = 0.98): a saturating fit explains 81% of the variance but underestimates the deepest cells, where density continues to rise to 170–189 kg m⁻³. Density also rises systematically with snow persistence, from 100–115 kg m⁻³ in cells with snow in fewer than 10% of months to 150–189 kg m⁻³ where snow is present in 45–65% of months, showing that the time-averaged density reflects the length of the snow season as well as compaction. The effective density implied by the SWE–SD regression slope (250 kg m⁻³) exceeds the median time-averaged density (123 kg m⁻³) and the median snow-month density (156 kg m⁻³) because long-term SWE is dominated by deep, dense snow. The correlations involving density depend strongly on the sampling unit: the SD–ρs coefficient falls from 0.87 (monthly spatial means) and 0.81 (cell means) to 0.62 (all cell-months), and within snow-present months its Spearman coefficient is only 0.58. This behaviour is consistent with the scale-dependent depth–density relationship reported by McCreight and Small (2014) and suggests that part of the inconsistency in the literature arises from differences in scale and sampling unit.
PC1 explains 90.8% of the variance, PC2 9.0% and PC3 0.2%. All variables load positively on PC1 (0.55–0.60), which therefore represents a general snow-accumulation gradient, whereas PC2 contrasts density (0.83) with depth and SWE. The elbow criterion indicated K = 4, the silhouette coefficient was highest at K = 2 (0.91; 0.74 at K = 4) and the Calinski–Harabasz index increased monotonically. K = 4 was retained because it separates the physically distinct mountain snowpack and is stable in the two-variable test (ARI = 0.91; 96.5% agreement). The classes form a continuum rather than discrete regimes: class 1 (snow-free to ephemeral; 59.6% of the area) holds 0.9% of the storage, class 2 (shallow seasonal; 32.1%) 25.7%, class 3 (persistent seasonal; 7.2%) 47.5% and class 4 (mountain snowpack; 184 cells, 1.2%) 25.8%.
The seasonal cycle has a December–March core with a February peak (median spatial-mean depth 0.062 m); its timing has remained stable, whereas its amplitude has roughly halved. December–March mean depth and SWE declined significantly (Z = −3.71 and −3.73; Sen's slopes −8.3 mm and −1.9 mm w.e. per decade), and the Pettitt test identified 2008 as the last year before a downward shift (p ≤ 0.002), after which the winter means fell by about half (SD 0.052 → 0.027 m; SWE 11.1 → 5.6 mm). April depth and SWE also declined significantly (change point 1998), whereas the November series showed no trend; the evidence therefore points to an earlier end, but not to a later start, of the snow season. Between the two periods, SWE decreased in 93.3% of the snow-affected area (median −5.1 mm; −41%). Cell-wise tests show significant winter decreases over 80.3% of that area and no significant increase; the largest losses occur on the flanks of the central Zagros, in the northwestern highlands and on the northern slopes of the Alborz, and March shows the most widespread significant decline (73.8% of active cells).
Conclusion
Within ERA5-Land, Iran's snow storage is concentrated in a few high-elevation cells and in a December–March window, and it has declined markedly since the late 2000s. Depth explains nearly all of the spatial variability of SWE, partly by construction, whereas density shows a narrow, nonlinear and sample-dependent relationship with depth that also encodes snow persistence. Two practical implications follow, both of which require in-situ validation. First, time-averaged reanalysis densities should not be used directly to convert depth to SWE; the effective density of 150–250 kg m⁻³ obtained here is more appropriate, and density models that vary in space and time are needed. Second, the concentration of storage in a small number of high-elevation cells, together with its recent decline, makes downstream water resources (Karun, Karkheh, Dez, Zayandeh-Rud and Sefidrud) vulnerable to further warming. The main limitations are the ~9 km resolution, the 100 kg m⁻³ density floor, the absence of independent validation, and spatial autocorrelation and multiple testing in the cell-wise maps.
Keywords
Subjects
