Large-ensemble climate simulations to assess changes in snow stability over northern Japan

Abstract

To examine the influence of global warming including increased heavy snowfall frequency on the potential of natural dry snow avalanche frequency and magnitude, we estimated the frequency of weak layer formation and the associated slab overload above the weak layer over northern Japan. The estimation was numerically performed using climate models' output for 1800 winter simulations in each of the historical (1951–2010) and +4°C experiments by forcing a physical-based snowpack model with the result of the climate models. Here the +4°C experiment was defined as a climate when the global mean air temperature had increased by 4°C from the preindustrial level. The estimation results showed that the probability of weak layer formation, identified by the natural stability index, would decrease all over the area because of the shorter age of the weak layers caused by a warmer climate, indirectly indicating a potential decrease in avalanche frequency. However, because of increased heavy snowfall frequency, slab overload would increase by 10–15% in inland areas for weak layers of decomposing fragments/precipitation particles and the mountainous area facing the Sea of Japan for weak layers of facets/depth hoar, thereby potentially indicating an increased magnitude of avalanches.

1. Introduction

Snow avalanches frequently destroy forests, infrastructures, houses and lives. Recent avalanche records showed a decreasing trend of avalanche occurrence, which was negatively/positively correlated with air temperature/snowfall, from the 1970s in the French Alps (Eckert and others, 2013) and Swiss Alps (Teich and others, 2012a) and from the 1950s in the US Northern Rocky Mountains (Peitzsch and others, 2021), especially in low elevations (Lavigne and others, 2015; Giacona and others, 2021). Similarly, the avalanche magnitude retrieved from runout distance or extent of avalanche-prone terrain has been decreased in recent decades below an elevation of 1200 m in the Vosges Mountains in France (Giacona and others, 2021). Nevertheless, an increase in wet avalanches would be expected because of a wetter snowpack resulting from a warmer climate (Hock and others, 2019). In particular, in high elevations, avalanches occurred frequently in recent decades in French Alps (Castebrunet and others, 2012). The frequency, and possibly the runout distance, and extent of avalanches retrieved from tree rings were also getting higher, longer and larger above an elevation of 2600 m in the Western Himalaya after the 1970s (Ballesteros-Cánovas and others, 2018), with a clear link to warming for the frequency increase. These recent changes were probably because snowpack would be getting wetter in midwinter and early spring, although the amount of snowpack in high elevation would be less sensitive to air temperature (López-Moreno and others, 2012; Sun and others, 2016; Katsuyama and others, 2017), or because snowfall intensity would be getting stronger mainly because of a moisture-rich atmosphere (Kawase and others, 2016; Sasai and others, 2019). Even at an elevation of 1200 m in the Vosges Mountains in France, avalanche activity remains in recent decades (Giacona and others, 2021). Hence, the spatial variation of avalanche formation and magnitude might be influenced by global warming, which potentially affects political/civilian strategies to adapt to natural avalanche hazards (Bebi and others, 2009; Strapazzon and others, 2021). To mitigate snow avalanche hazards, it is necessary to assess the effect of global warming on natural avalanches.

Martin and others (2001), to the best of our knowledge, were the first to attempt to assess the potential impact of global warming on avalanches using a stability index (SI) calculated using a physical snowpack model imposed with perturbed atmospheric data. Then, they found that owing to the warming climate, the potential occurrences of avalanches in the French Alps could decrease in winter and increase in spring. Notably, the estimation was based on the low stability of snowpack diagnosed by an SI, traditionally used to diagnose snowpack mechanical stability (Föhn, 1987; Jamieson and Johnston, 1998), which is a necessary condition for avalanche release (Reuter and Schweizer, 2018). Richter and others (2020) also performed a similar numerical experiment to that by Martin and others (2001) and concluded that the snowpack stability diagnosed by the SI increased and decreased with the increase in air temperature and precipitation, respectively. Further, Castebrunet and others (2014) estimated the global warming impact on avalanche activity by establishing a statistical relationship between avalanche activity and snowpack instability in the meteorological state of the present climate and by seeking a similar meteorological state between the present and future climates. Then, they determined that avalanche activities decreased and increased at low and high altitudes, respectively. However, their statistical approach requires enormous past avalanche records. Because there is no such dataset in many countries, including Japan, the numerical model approach using the SI is the only feasible technique to assess the impact of global warming at present.

Notably, in the assessment, increases in air temperature and heavy snowfall could have contrasting effects on snowpack mechanical stability. High air temperatures could improve the shear strength of snow layers in a shorter period than low air temperatures due to faster metamorphism (Lehning and others, 2002), which could increase stability and shorten the lifetime of weak layers. This effect of temperature increase would also reduce the amount of overlying slab above the weak layer because the slab was fed by less amount of snowfall during the shorter lifetime of the weak layer. However, an increase in heavy snowfall frequency could rapidly increase the amount of the overlying slab and could decrease snowpack stability. Moreover, a thin snow depth due to a warmer climate would serve as a preferable environment for weak layer formation because of the high vertical gradient of snow temperature (Birkeland and others, 1998). To elucidate these complex effects of global warming on snowpack stability, researchers have employed physical-based snowpack models, such as SNOWPACK (Bartelt and Lehning, 2002) and Crocus (Brun and others, 1992), which produce a time evolution of a vertical snowpack layer structure with snowpack parameters related to avalanche formation, such as hardness, density, grain type, shear strength and shear stress.

From the perspective of the effects of increased heavy snowfall frequency on avalanche occurrence, numerical snowpack simulations with large ensembles of climate data could be used because the frequency of extreme events is usually low. Recently, Mizuta and others (2017) performed large-ensemble climate simulations and arranged their results as a dataset named the database for policy decision making for future climate change (d4PDF). This dataset stores the results of large-ensemble simulations dynamically downscaled from the results of a general circulation model (GCM) with a regional climate model (RCM) with a 20 km spatial resolution for 3000 years of historical climate and 5400 years of future climate in total. By analyzing d4PDF, Kawase and others (2016) found that the frequency of heavy snowfall would increase in mountainous areas in central Japan and would decrease in coastal areas along the Sea of Japan. However, the avalanche affected by the increase in heavy snowfall has not been addressed.

In this study, we identified the impact of global warming on the potential of natural dry snow avalanches over northern Japan, including the effect of increased heavy snowfall, by forcing the SNOWPACK model with d4PDF. The procedure of the SNOWPACK model calculation was nearly identical to that used by Katsuyama and others (2020), who estimated the effect of global warming on the amount and grain types of snowpack over Hokkaido, Japan's northern island, by forcing the SNOWPACK model with future climate variables simulated by an RCM with a 10 km spatial resolution. Avalanche potential was indicated by the SI, defined as a ratio of shear strength and shear stress. Even though the approach using the SI could only diagnose a necessary condition for avalanches (i.e. not direct estimation of the global warming effect on avalanches), this would be the only feasible technique because of no dataset containing enormous past avalanches in Japan. Because of the same reason, we did not address the physical/statistical relation between the avalanche potential and the actual avalanche activity as it should be addressed in the future work when enough avalanche records are gathered. Moreover, we did not address wet avalanches because the process of wet avalanche initiation could be related to the interactions between snow grains and liquid water in the snowpack rather than the mechanical characteristics, such as shear strength and shear stress (e.g. Mitterer and others, 2011).

The remainder of this article is organized as follows. In Section 2, the basic atmospheric and snowpack characteristics of the study area are presented. The long-term observation data and d4PDF are briefly described in Section 3. The methodology used to estimate avalanche potential over northern Japan using d4PDF and the SNOWPACK model is described in Section 4. In Section 5, we first briefly check the changes in snowfall including their extremes given by d4PDF. Then, the snow water equivalent (SWE), snow-covered days (SCD) and snow grain types calculated using the SNOWPACK model in the present and future climates are presented to check whether our results agree with those of previous studies. Further, the long-term variability of snowpack amount determined via the calculation is compared with that determined via observation. Then, changes in snowpack properties potentially related to avalanche occurrence and magnitude between the present and future climates are shown. In Section 6, we first check the consistency of the obtained results with associated previous studies. Then, we provide interpretation and limitations of our results and the challenges that future research must address. Finally, in Section 7, we conclude this study.

2. Study area

In this study, we focused on northern Japan, which features southwest–northeast mountains in Honshu and Hokkaido islands (Fig. 1). In winter, there is a considerable amount of snowfall in the coastal region facing the Sea of Japan because winter monsoonal winds from northwest to southeast advect much of the vapor produced by the Sea of Japan (e.g. Takano and others, 2008). In particular, in the Hokuriku region (Fig. 1), snowfall is heavier than that in other regions because of an updraft enhanced by the mountains (Houze, 2014) and moisture flux convergence in the central Sea of Japan, usually called the Japan Sea Polar Air Mass Convergence Zone (Nagata, 1987), whereas the air temperature typically holds steady at 0°C, except in the mountains. In central Honshu, which is located leeward of the winter monsoonal winds, snowfall occurs primarily because of extratropical cyclone passage rather than winter monsoonal winds. Because of this different snowfall mechanism, central Honshu has less snowfall than Hokuriku. In Hokkaido, the northern island of Japan (Fig. 1), which is characterized by a low air temperature, there is a large snowfall on the western side due to winter monsoonal winds and less snowfall on the eastern side due to extratropical cyclone passage.

Fig. 1. (a) Topography of the Meteorological Research Institute atmospheric general circulation model (MRI-AGCM) and the nonhydrostatic regional climate model (NHRCM). The solid black and red lines indicate the calculation domains of the NHRCM and SNOWPACK, respectively. (b) Locations of Asahikawa, Tokamachi and Kiso (red type) and regions of Hokkaido, Hokuriku and central Honshu (blue outlines) depicted with the topography of the NHRCM.

Asahikawa, Tokamachi and Kiso were selected as representative locations of central Hokkaido, Hokuriku and central Honshu, respectively (Fig. 1b). Asahikawa is characterized as having a low air temperature and moderate snowfall amount in the study area. During 1991–2020, the climatological precipitation, air temperature and maximum snow depth in December–February (DJF) were 222.4 mm, −5.7°C and 85 cm, respectively (JMA, 2021). Because of this climate, the snowpack mainly comprises dry snow and sometimes faceted crystals and depth hoar (FC/DH) (Ishizaka, 2008; Katsuyama and others, 2020). Tokamachi is characterized as having a large amount of snowfall, a typical air temperature of 0°C, and a relatively large number of SCD. During 1918–2017, the climatological precipitation, air temperature, maximum snow depth and the number of SCD in DJF were 1053 mm, 0.8°C, 230 cm and 116 d, respectively (Takeuchi and others, 2019a). The relatively high air temperature rapidly metamorphoses precipitation particles and decomposing fragmented precipitation particles (PP/DF) to wet snow, whereas the frequent snowfall often results in the formation of PP/DF snow layers at the top of snowpacks (Takeuchi and others, 2019b). Kiso is characterized as having a slightly low air temperature and less snowfall. During 1991–2020, the climatological precipitation and air temperature in DJF were 220.7 mm and −0.2°C, respectively (JMA, 2021). Hence, FC/DH is primarily formed in Kiso (Ishizaka, 2008).

3. Data

3.1. Long-term snowpack observation

Long-term observation of SWE at Tokamachi, Niigata, Japan (Fig. 1) (Takeuchi and others, 2019a), was used to confirm the probability density function (PDF) of the seasonal maximum SWE (SM-SWE) in the present climate calculated using the SNOWPACK model forced with d4PDF. The observation of SWE was performed over 81 years, beginning in the 1939/1940 winter. Data were manually collected every several days between 9:00 and 10:00 Japanese local time, but the interval depended on the periods (Table 1). Starting in 2005, automatic SWE measurements taken using a snow weight meter at 9:00 were used instead of data from manual observation. The PDF of SM-SWE was constructed from a series of observed SM-SWEs.

Table 1. Observation intervals of SWE for each period at Tokamachi, Niigata, Japan

3.2. d4PDF

In this study, the SNOWPACK model was forced with the results of the RCM in d4PDF (Mizuta and others, 2017). These data were created by dynamically downscaling using a nonhydrostatic RCM (NHRCM) with a 20 km spatial resolution from version 3.2 of the Meteorological Research Institute atmospheric GCM with a 60 km spatial resolution. The calculation domain of the NHRCM includes all of Japan's islands (Fig. 1a). In the historical climate experiment, 50 ensemble simulations were performed from 1 September 1950, to 31 August 2011, using different initial conditions and the observed perturbated sea surface temperatures (SSTs). Here, the warming level of the global mean air temperature during 1950–2011 is +0.45°C from the preindustrial level (Mizuta and others, 2017). The future climate simulation was performed under conditions in which the global mean air temperature increases by 4°C from the preindustrial condition under the Representative Concentration Pathway 8.5 scenario (i.e. the increase in the global mean air temperature from the historical level was +3.55°C). In this +4°C experiment, 15 ensemble simulations were performed using different initial conditions and perturbated SST changes projected by six different GCMs of the Coupled Model Intercomparison Project Phase 5. In this study, owing to the limited computational resources, we only selected data from 30 ensemble members collected from 1 August 1951 to 31 July 2011, in the historical experiment and data from five ensemble members with six different GCMs collected from 1 August 2051 to 31 July 2111, in the +4°C experiment. The total data length was 1800 years in each of the historical and +4°C experiments (Table 2). Notably, the global warming level was adjusted to keep the homogeneous level of +4°C throughout the period from 2051 to 2111 (Mizuta and others, 2017). Moreover, it should be noted that, in the ensemble experiment like this study, the hourly (or daily) time series are not comparable to the observations, but are comparable in terms of distribution (or associated statistics such as mean and std dev.).

Table 2. Selected d4PDF

^a Community Climate System Model version 4 (CCSM4).

^b Geophysical Fluid Dynamics Laboratory Climate Model version 3(GFDL CM3).

^c Hadley Centre Global Environment Model version 2 (HadGEM2).

^d Model for Interdisciplinary Research on Climate version 5 (MIROC5).

^e Max Planck Institute Earth System Model (MPI-ESM).

^f Meteorological Research Institute Coupled Atmosphere-Ocean General Circulation Model version 3 (MPI-CGCM3).

The NHRCM has been confirmed to reproduce the spatial snowfall pattern and basic weather pattern that result in heavy snowfall over Japan, as introduced in Section 2 (Kawase and others, 2016), including snowfall that occurs when the air temperature is ~0°C, which is a characteristic of Hokuriku (Ohba and Sugimoto, 2020). However, compared with automatic weather stations maintained by the Japan Meteorological Agency, the frequency of heavy snowfall by the NHRCM was underestimated in the low altitude of the coastal region along the Sea of Japan and overestimated in the midaltitude mountainous region along the Sea of Japan (Kawase and others, 2016).

4. Method

4.1. SNOWPACK model

We performed snowpack calculations by forcing the results of the NHRCM to the SNOWPACK model (version 3.6.0), which can provide complete hourly time series of snowpacks' structures and physical properties, such as grain type, density and shear strength (Bartelt and Lehning, 2002), for all ensemble members of both the historical and +4°C experiments (Supplement 1 provides some example time series of SWE at Tokamachi). The calculation period for the historical experiment was from 1951/1952 winter to 2010/2011 winter, and the calculation period for the +4°C experiment was from 2051/2052 winter to 2110/2111 winter. The calculations were started on August 1 of each year and stopped on July 31 of the following year. The SNOWPACK model was forced with seven atmospheric variables (air temperature, wind speed, relative humidity, incoming shortwave and longwave radiations, and solid and liquid precipitations) every hour. The model was initialized with 10 m of soil layers with a temperature of $0.97\bar{T} + 2.3\,[ ^\circ {\rm C] }$°C, where $\bar{T}$ is the climatology of the annual mean air temperature (Hirota and others, 1995), and the soil temperature of the bottom boundary was fixed at the initial temperature. Although SNOWPACK was originally a 1-D model for a specific site, we calculated the spatial distribution of snowpacks simply by looping the model at each gridpoint. The grid configuration and elevation at each gridpoint were the same as those of the NHRCM, but the calculation domain was smaller to cover only northern Japan (Fig. 1a). The model calculation was performed by assuming that each gridpoint represented flat terrain without vegetation so that the effects of snow redistribution by wind and the modification of meteorology by vegetation would not be included. The possible problem related to this assumption will be discussed in Section 6. See Katsuyama and others (2020) for details regarding the configuration of the SNOWPACK model.

When assessing the impact of global warming, it is imperative to check the validity of the historical results represented by the numerical model. Fortunately, the SNOWPACK model has been well confirmed to reproduce reasonable snowpack structures, including their physical properties, by comparing the calculation results and snow pit observations from many countries, including Japan (e.g. Bartelt and Lehning, 2002; Hirashima and others, 2004, 2008, 2009; Yamaguchi and others, 2004; Nishimura and others, 2005; Katsuyama and others, 2020). Hence, we did not validate the simulation skills of SNOWPACK itself in this study. However, because this is the first attempt to apply large-ensemble data simulated by climate models to the SNOWPACK model in Japan, the PDF of SM-SWE of large-ensemble simulations was compared with the long-term SWE observation at Tokamachi (Subsection 3.1). In this comparison, we selected three model grid points nearest to Tokamachi because it is located at an intermediate point among the grid points. Then, the PDF constructed from all seasonal series of the calculated SM-SWE at the three grid points was compared with the PDF constructed from the seasonal series of the observed SM-SWE. Notably, the calculations based on climate models are comparable to the observations only in terms of distribution (or associated statistics such as mean and std dev.). Therefore, for example, the daily time series of snow depth in a specific year were not correlated between the calculation and observation at all (Supplement 2).

4.2. Fundamental snowpack characteristics

In this study, we computed some simple statistics from the hourly time series of the SNOWPACK model outputs at each gridpoint: SM-SWE, the number of SCD in a season where snow depth was more than 0 cm, and a thickness fraction of grain type, a simple statistic that indicated the fraction of snow grain type throughout a season, defined as the fraction of the time-integrated layer thickness of a particular snow grain type to the time-integrated snow depth during a season (Katsuyama and others, 2020). These variables were computed grid-by-grid independently of each winter and ensemble member. The SWEs and SCD were averaged over each of the historical and +4°C experiments as an arithmetical mean. An average thickness fraction across the ensembles and winter is a weighted mean, $\bar{x} = \sum\nolimits_i {f_id_i/} \sum\nolimits_i {d_i}$, where f is the thickness fraction and d is the time-integrated snow depth during a season. A future change in the mean of SM-SWE was assessed through a relative difference between the values of the historical and +4°C experiments to that of the historical experiment. A future change in the mean of SCD was simply shown by an absolute difference between the values of experiments (future minus historical).

4.3. Diagnosis of avalanche potential

Dry snow avalanche potential is hourly calculated based on an index generated from the results of the SNOWPACK model. Many operational avalanche warning services diagnose the avalanche potential based on indices (Morin and others, 2020). A natural SI is considered the simplest index and is often used to estimate the potential of avalanches, assuming a uniform slope (Föhn, 1987; Jamieson and Johnston, 1998):

(1)$$S_\theta = \displaystyle{\sigma \over {w\cos \theta \sin \theta }}, \;$$

where S is the SI, w [Pa] is the vertical stress due to the weight of the overlying slab above a specific snow layer, θ [°] is the slope angle and σ [Pa] is the shear strength of the snow layer. A low index value indicates that the snow layer is too mechanically weak to support the overlying snow layers. Thus, the snow layer with the minimum SI is mechanically the most frangible (i.e. the weakest), and the snow layers above the weak layer are released by a failure. In this study, SI was calculated by assuming a 38° slope, S ₃₈, in Eqn (1) because that is the most common avalanche slope angle (McClung and Schaerer, 2006). Although several other indices can be used to estimate the potential of avalanches, such as the skier SI (Föhn, 1987; Jamieson and Johnston, 1998; Schweizer and others, 2006; Monti and others, 2016) and the structural SI (Schweizer and others, 2006), they can only calculate the potential of human-triggered avalanches, not natural avalanches.

The weak layer was time independently defined with a distinction of two different types depending on the snow grain types of PP/DF and FC/DH. These grain types often act as a weak layer through the different physical processes: PP/DF becomes a weak layer when snowfall intensity is strong relative to strengthening the speed of snow grains by sintering (Endo, 1993), and FC/DH becomes a weak layer when the shear strength is weakened by kinetic growth metamorphism (Hirashima and others, 2009). The PP/DF and FC/DH weak layers were, respectively, identified on the basis of where on the PP/DF and FC/DH snow layers the minimum SI was <1.5, which was considered the threshold value because many natural avalanches occurred when the SI was <~1.5 and this value has traditionally been used to identify unstable snowpacks (Perla, 1977; Jamieson and Johnston, 1995; Lehning and others, 2004; Hirashima and others, 2008). This threshold method using SI estimates the best probability of capturing a necessary condition for natural dry snow avalanche formation (i.e. the best probability of detecting nonavalanche days) in Switzerland (Reuter and others, 2022). Notably, two weak layers would be defined in maximum for a single time step when both PP/DF and FC/DH are included in snowpacks because the weak layers of PP/DF and FC/DH are independent. Finally, the potential magnitude of avalanches was identified on the basis of the amount of the slab overlying the weak layer. The weak layer and slab overload, as determined using SI, were analyzed grid-by-grid, in elevation ranges, and at the three locations (Asahikawa, Tokamachi and Kiso) (Fig. 1b). Moreover, according to Gaume and others (2015), when the slab density exceeds 150 kg m⁻³, the potential release zone for dry slab avalanches would theoretically extend to the potential maximum zone mainly controlled by the topography. Since the potential magnitude of avalanches should be proportional to the multiplies of the amount of slab and the area of the release zone, we also addressed the slab density in addition to the amount of slab at the three locations.

In the grid-by-grid and elevation range analyses, the identified weak layers were further processed into the probability normalized with SCD in each of the historical and +4°C experiments:

(2)$$P \vert _{S_{38} < 1.5} = \mathop \sum \limits_{i = 1}^{N_m} \displaystyle{{\mathop \sum \nolimits_{\,j = 1}^{24{N_{s}}_{i} } H_0( {1.5-{S_{38}}_{i, j} } ) } \over {24{N_{s}}_{i} }}, \;$$

where N_m (=1800) is the total number of winters of each experiment, N _s is the number of SCD and H ₀ is the Heaviside step function. The probability at each elevation range is the mean weighted by the SCD at the gridpoint corresponding to the elevation range. The slab overloads corresponding to the identified weak layers were also processed into the mean value for all weak layers for all ensemble members in each of the historical and +4°C experiments. In the analysis for the three locations, we selected the three nearest grid points on the model to a location. Then, the PDFs of SI, slab overload and slab density at the locations were constructed using all hourly time series of each of the historical and +4°C experiments, but nonsnow-covered times were excluded. Similarly, a joint PDF of slab overload and weak layer age from the layer birth was also constructed. The model elevations of the selected grid points were 231, 309 and 280 m in Asahikawa; 347, 537 and 767 m in Tokamachi; and 1330, 1297 and 1199 m in Kiso (Fig. 1).

5. Results

5.1. Changes in snowfall

First, we check the future changes in snowfall produced by d4PDF including heavy snowfall events. Figure 2 shows the cumulative frequency distribution (CFD) of daily snowfall (without liquid precipitation) from November to the following March (NDJFM) at Asahikawa, Tokamachi and Kiso. In Asahikawa, between the historical and +4°C experiments, the probabilistic snowfall of once per winter, once per 10 winters and once per 100 winters all increased slightly by ~2 kg m⁻² d⁻¹ (Fig. 2a). This pattern of change in CFD directly means the increase in the frequency of heavy snowfall. However, even with the increases in the heavy snowfall, the mean daily precipitation in NDJFM would not change (figure not shown). In Tokamachi, the snowfall amount with a longer/shorter return period than once per 5 winters, that was ~75 kg m⁻² d⁻¹, would increase/decrease in the +4°C experiment (Fig. 2b). Even if the frequency of heavy snowfall would increase, however, the total snowfall amount in NDJFM in the +4°C experiment would decrease by ~45% comparing that in the historical experiment (figure not shown) because of warmer air temperature. In Kiso, the probabilistic snowfall amount of once per winter, once per 10 winters and once per 100 winters all decreased slightly (Fig. 2c), meaning the increase in the frequency of weak snowfall intensity. Hence, the total snowfall amount in NDJFM would decrease by ~25% in the +4°C experiment (figure not shown).

Fig. 2. Cumulative frequency distribution of daily snowfall amount from November to March in the following year at (a) Asahikawa, (b) Tokamachi and (c) Kiso. The black and red lines indicate the result of the historical and +4°C experiment, respectively. The return period of once per winter, once per 10 winters and once per 100 winters is denoted by a dotted line and a text, respectively.

5.2. SWE, SCD and grain types

Figure 3 shows the mean of the SM-SWE in the historical experiment and its difference from that of the +4°C experiment. The spatial distribution of the mean value was characterized as having a high SWE along the Sea of Japan and in central Hokkaido (Fig. 3a). The SWE was especially high in the mountainous area of Hokuriku because of the heavy snowfall (Section 2; Fig. 2b). In Hokkaido, the SWE was ~600 kg m⁻² in the coastal region along the Sea of Japan in the historical experiment (Fig. 3a). However, in central Hokkaido, the SWE was ~1000 kg m⁻² in the historical experiment. Global warming would largely decrease the mean of the SM-SWE over northern Japan (Fig. 3b). In several areas, except the mountainous area of Hokuriku and most areas of Hokkaido, the SWE would decrease by ~80% compared with the historical experiment. In mountainous areas, the decrease would be relatively small, i.e. ~20–30%.

Fig. 3. (a) Mean of the seasonal maximum snow water equivalent (SM-SWE) [kg m⁻²] of the historical experiment, (b) relative difference of the mean SM-SWE between the +4°C and historical experiments normalized by the SWE of the historical experiment and (c) probability density function (PDF) of SM-SWE at Tokamachi. In (c), the gray bars, solid black line and dashed black line indicate the observation, the result of the historical experiment and the result of the +4°C experiment, respectively.

The PDF of the SM-SWE was also compared with that of Tokamachi (Fig. 3c). It showed a broad distribution with a mean of 800 kg m⁻² and a range of 100–1500 kg m⁻². The PDF reproduced by the SNOWPACK model forced with d4PDF in the historical experiment was similar to the observed PDF. Global warming largely decreased the SWE at all ranges of values: the most frequent value in the +4°C experiment was ~100 kg m⁻², corresponding to the minimum value of the historical observation (Fig. 3c). Even the top-10th percentile value in the +4°C experiment (560 kg m⁻²) was less than the mean value in the historical experiment. Notably, the double-peak shape of the PDF in the historical experiment was simply because multiple grid points on the model were combined (Subsection 4.1).

Figure 4 shows SCD during a season. In the historical experiment, in Hokkaido, there were 180 SCD in the central mountainous region and ~150 SCD in other regions (Fig. 4a). In Hokuriku, there were more than 180 SCD in the mountainous region, and the number decreased closer to the Sea of Japan. In central Honshu, there were ~60–150 SCD. Global warming would monotonically decrease the number of SCD over northern Japan (Fig. 4b). In Hokkaido, the number of SCD would uniformly decrease by 30 d. In the mountainous area of Hokuriku and central Honshu, the SCD would decrease by ~60 d. In the coastal area of Hokuriku, the decrease would be more than 80 d, which mostly indicated no snow-covered.

Fig. 4. (a) Mean of seasonal snow-covered days (SCD) of the historical experiment and (b) difference in SCD between the +4°C and historical experiments.

Figure 5 shows the thickness fraction of grain types. The fraction of PP/DF largely increased from ~20% in the historical experiment to ~40% in the +4°C experiment, except in Hokkaido (Figs 5a, c). In Hokkaido, the fraction of PP/DF would increase from ~10 to ~20%. FC/DH was estimated to exist mainly in eastern Hokkaido and slightly in western Hokkaido and central Honshu in the historical experiment (Fig. 5b). In the +4°C experiment, FC/DH would decrease, especially in Hokkaido (Fig. 5d).

Fig. 5. (a) Thickness fraction of precipitation particles and decomposing and fragmented precipitation particles (PP/DF) and (b) that of faceted crystals and depth hoar (FC/DH) relative to the time-integrated snow depth during a season for the historical experiment. (c) and (d) are, respectively, identical to (a) and (b), except for the +4°C experiment. Areas, where the SM-SWE was <50 kg m⁻² in the historical experiment, were masked out (Fig. 3a).

5.3. Weak layer formation

Figure 6 shows the existence probabilities of weak layers where the minimum SI value across the snowpack is <1.5 (Subsection 4.3) during the snow-covered duration. In the historical experiment, weak layers of PP/DF frequently formed in several areas corresponding to the spatial distribution of SWE (Fig. 3a): the probability of weak layer formation was more than 20% in the mountainous region of Hokuriku and ~10% in Hokkaido. Following this spatial pattern, the probability was well correlated with elevation (Fig. 7a). However, weak layers of FC/DH exist mainly in Hokkaido and central Honshu (Fig. 6b) in correlation with the FC/DH fraction (Fig. 5b) rather than the SWE (Fig. 3a). Hence, the probability did not correlate with elevation (Fig. 7b).

Fig. 6. Probability of weak layer formation identified by the threshold of the natural stability index (SI) for (a) PP/DF and (b) FC/DH in the historical experiment and (c) PP/DF and (d) FC/DH in the +4°C experiment. Areas, where the SM-SWE was <50 kg m⁻² in the historical experiment, were masked out (Fig. 3a).

Fig. 7. Probability of (a) PP/DF and (b) FC/DH weak layers at each elevation range in the (blue bar) historical and (red bar) +4°C experiments. Areas, where the SM-SWE was <50 kg m⁻² in the historical experiment, were excluded (Fig. 3a). The bins for elevation are 100 m intervals.

Moreover, because the historical experiment includes results regarding the snowpack calculation for 1800 winter simulations (Subsection 3.2), we can address the PDF of the minimum SI value across the snowpack, which may represent the climatological characteristics of snowpack stability at a specific site. The PDF of SI shows a lognormal-like distribution rather than a normal distribution independent of locations and grain types (Fig. 8). In Asahikawa in the historical experiment, the SI value for both PP/DF and FC/DH was the most frequent at 2, which was extremely close to the threshold of the weak layer (Subsection 4.3). In Tokamachi, the most frequent SI value was 1.5, a threshold value of the weak layer, for both PP/DF and FC/DH. In Kiso, the most frequent SI value was ~3.5 for both PP/DF and FC/DH. Notably, the PDF of the SI value integrated between 0 and 1.5 is not equal to the probability of weak layer formation normalized with the SCD during a season (Fig. 6) because PP/DF or FC/DH is not always included in the snowpack for every time step.

Fig. 8. The PDF of SI in Asahikawa (red), Tokamachi (green) and Kiso (blue) in the historical (solid line) and +4°C (dashed line) experiments for (a) PP/DF and (b) FC/DH.

Global warming would decrease the probability of weak layer formation for both PP/DF and FC/DH (Fig. 6) at all elevation ranges (Fig. 7). The probability for PP/DF was large only in the mountainous area of Hokuriku and central Hokkaido in the +4°C experiment (Fig. 6c). For FC/DH, in all regions, the weak layer would rarely form (Fig. 6d). The peak of the PDF of the SI shifted to a stable side, typically for PP/DF and FC/DH in the three locations (Fig. 8).

5.4. Slab overloads and density

The slab overloads above the weak layer that coincide with the SI threshold of 1.5 were further addressed by the spatial distribution of the mean value (Fig. 9) and its dependency on elevation (Fig. 10), and in the three locations, the PDFs of the slab overload and density were also addressed (Figs 11–13). In the historical experiment, the mean overload above the weak layer of PP/DF was ~100 kg m⁻² in the mountainous region of Hokuriku, which was considerably larger than that of the other regions, ~40 kg m⁻² in Hokkaido and the central Honshu (Fig. 9a). The spatial distribution of overload for FC/DH was similar to that for PP/DF in the historical experiment because the overload was correlated with snowfall intensity (Fig. 9b). However, because the weak layer of FC/DH would last a relatively long time depending on the temperature gradient of the snowpack and on how much snowfall accumulates on the weak layer, the amount of overload for FC/DH was much larger than that for PP/DF (Figs 9a, b). In the mountainous region of Hokuriku, in particular, where there was a considerably large amount of snowfall, the amount of overload for FC/DH would be ~160 kg m⁻², much larger than that for PP/DF. However, in Hokkaido and central Honshu, the overload for FC/DH in the historical experiment was ~40 kg m⁻², mostly the same amount as that for PP/DF (Figs 9a, b). Following this spatial pattern, the overloads for both PP/DF and FC/DH were well correlated with elevation (Fig. 10).

Fig. 9. Slab overload above the weak layer identified with an SI threshold <1.5 for (a) PP/DF and (b) FC/DH in the historical experiment and its relative difference between the historical and +4°C experiments for (c) PP/DF and (d) FC/DH normalized by the slab overload of the historical experiment. Areas, where the SM-SWE was <50 kg m⁻² in the historical experiment, were masked out (Fig. 3a).

Fig. 10. Slab overload for (a) PP/DF and (b) FC/DH weak layers at each elevation range in the (blue bar) historical and (red bar) +4°C experiments. Areas, where the seasonal maximum SWE was <50 kg m⁻² in the historical experiment, were excluded (Fig. 3a). The bins for elevation are 100 m intervals.

Fig. 11. PDF of slab overload for (a, c, e) PP/DF and (b, d, f) FC/DH at (a, b) Asahikawa, (c, d) Tokamachi and (e, f) Kiso. The black and red lines indicate the results of the historical and +4°C experiments, respectively.

Fig. 12. Joint PDF of slab overload and the weak layer age from the layer formation in (a, b) Asahikawa, (c, d) Tokamachi and (e, f) Kiso, respectively. (a, c, e) and (b, d, f) indicate the weak layer of PP/DF and FC/DH, respectively. The black and red contour lines indicate the historical and +4°C experiments, respectively. Thick and thin contours show probability density values of 10⁻³ and 10⁻⁴, respectively.

Fig. 13. PDF of slab density for (a, c, e) PP/DF and (b, d, f) FC/DH at (a, b) Asahikawa, (c, d) Tokamachi and (e, f) Kiso. The black and red lines indicate the result of the historical and +4°C experiment, respectively. The black and red text beside the gray vertical dashed line shows the probability integrated over the density larger than 150 kg m⁻³ in the historical and +4°C experiment, respectively.

Figure 11 shows the PDF of the slab overload in the three locations. The overload for PP/DF in the historical experiment was the most frequent at ~50 kg m⁻² in Asahikawa (Fig. 11a) and Tokamachi (Fig. 11c) and at ~20 kg m⁻² in Kiso (Fig. 11e). The top-10th percentile value was approximately two or three times larger than the most frequent value: 85, 120 and 60 kg m⁻² in Asahikawa, Tokamachi and Kiso, respectively. Similarly, the projected value range of overload for FC/DH is broad in all three locations in the historical experiment (Figs 11b, d, f). Specifically, in Tokamachi, the overload would range from ~0 to 300 kg m⁻², the most frequent value being 150 kg m⁻² (Fig. 11d).

Because of global warming, the mean slab overload would increase in central Honshu and slightly increase in central Hokkaido for the PP/DF weak layer (Fig. 9c) and slightly increase in the mountainous region of Hokuriku and central Hokkaido for the FC/DH weak layer (Fig. 9d), whereas that would decrease in many other areas (Fig. 9). In central Honshu, in particular, the mean slab overload for PP/DF would increase by ~30% relative to the historical experiment (Fig. 9c). Because the PDF peak of slab overload for PP/DF would significantly shift from 10 to 50 kg m⁻² in Kiso (Fig. 11e), the increase in the mean slab overload in central Honshu would occur simply because of the peak shift. In Asahikawa, the probability of slab overload for 20 kg m⁻² in the historical experiment would be half in the +4°C experiments (Fig. 11a), and this would cause a slight increase in mean slab overload in central Hokkaido (Fig. 9c). Because these increases in the overload for PP/DF would be estimated in some limited areas, in terms of the dependency on elevation, no obvious change in the overload was seen (Fig. 10a). The mean overload for FC/DH would decrease in several areas, but a slight increase would be projected in Hokuriku and central Hokkaido. This increase in mean slab overload would be due to an increase in probability corresponding to a 180 kg m⁻² slab overload in Tokamachi (Fig. 11d) and that corresponding to a 100 kg m⁻² slab overload in Asahikawa (Fig. 11b) in the PDFs between the historical and +4°C experiments. In terms of the dependency on elevation, the overload for FC/DH would increase/decrease at the elevation range of 800–1500/0–500 m (Fig. 10b).

Because the slab overloads should be correlated with the snowfall amount from the snow layer initiation to the weak layer formation (Section 1), joint PDFs of the age of the weak layer and slab overloads were addressed (Fig. 12). In Asahikawa, the age mainly ranged in 0–10 d for the PP/DF weak layer (Fig. 12a) and 3–12 d for the FC/DH weak layers (Fig. 12b) in the historical experiment. In Tokamachi, the age ranges were 0–4 and 5–8 d for PP/DF and FC/DH, respectively. In Kiso, the ranges were 0–4 and 2–7 d for PP/DF and FC/DH, respectively. Commonly in the three locations, the slab overloads increased with an increase in weak layer age. Global warming would commonly reduce the age of the weak layers by a few days in the three locations for both weak layers (Fig. 12). In addition, the increasing speed of slab overload in the +4°C experiment would be faster than that in the historical experiment.

The changes in the slab overloads and density should be the two sides of the same coin via the settlement process. In Asahikawa, following the increase in the probability of the overload for 50 kg m⁻² for PP/DF weak layers due to the global warming (Fig. 11a), the most frequent density of slab in the +4°C experiment increased by ~10 kg m⁻³ compared to the historical experiment (Fig. 13a). Similarly, the peak of PDF of the slab density for FC/DH weak layers would shift to the larger density side (Fig. 13b). In Tokamachi, while the change in the PDF for PP/DF would be small between the historical and +4°C experiments (Fig. 13c), that for FC/DH would largely shift to the denser side (Fig. 13d). Following this, the probability integrated over 150 kg m⁻³, the probability which the avalanche release zone fully extended to the potential maximum zone mainly controlled by the topography (Gaume and others, 2015), was increased from 44% in the historical experiment to 76% in the +4°C experiment. In Kiso, following the increase in the probability of the overload for 50 kg m⁻² for the PP/DF (Fig. 11e), the PDF of slab density would slightly shift to the denser side (Fig. 13e) while that for FC/DH would not change mostly (Fig. 13f).

6. Discussion

6.1. Consistency of SWE, SCD and grain types

At the beginning of this study, we estimated the SM-SWE, SCD and thickness fraction of grain types in both the historical and +4°C experiments. Because some parts of these variables had already been estimated by previous studies, the consistency of our calculation could be confirmed simply by comparison with the previous results. The SWE spatial pattern in the historical experiment (Fig. 3a) was similar to that found by Hosaka and others (2005), who estimated the SWE over northern Japan using a 20 km RCM and a simple land surface model. In Hokkaido, Katsuyama and others (2020) estimated that the SWE was ~600 kg m⁻² in the coastal region along the Sea of Japan by a combination of a 10 km RCM and SNOWPACK model, which agreed with our results. However, in central Hokkaido, the SWE was ~1000 kg m⁻² in the historical experiment, which was ~150% heavier than the estimation using a 10 km RCM (Katsuyama and others, 2020). This was probably because the 20 km spatial resolution of the NHRCM was insufficient to represent the mountainous region, thereby overestimating the snowfall amount. For future changes, several previous studies reported that a decrease in SWE responding to global warming would be relatively small in mountainous areas (Hosaka and others, 2005; Katsuyama and others, 2020) because the amount of snowpack at a high altitude was correlated with precipitation rather than air temperature (Sun and others, 2016). The same decreasing pattern was obtained in this study. The monotonical decrease in SCD by ~30 d in Hokkaido agreed well with the result of Katsuyama and others (2020).

The thickness fraction of PP/DF would largely increase, except in Hokkaido, in this study because of the warmer climate (Figs 5a, c). This was because snow grains metamorphosed from PP/DF would disappear earlier than those in other regions because of the fewer SCD (Fig. 4b) due to the warmer climate (Katsuyama and others, 2020). The spatial pattern of FC/DH in the historical experiment (Fig. 5b) agreed well with the estimation of primary grain type over Japan based on a statistical relationship between atmospheric climatology and grain types observed (Inoue and Yokoyama, 2003; Ishizaka, 2008).

6.2. Weak layer formation and associated slab overloads and density

In this study, we estimated the decrease in weak layer formation over northern Japan (Figs 6, 8) and the increase in slab overload in some limited areas (Figs 9, 11) responding to global warming. These results could be interpreted by recalling the contrasting effect of the changes in air temperature and snowfall intensity on the avalanche potential: the increase in air temperature rapidly strengthens the snowpack and decreases the slab overload and the increase in precipitation intensity weakens the snowpack and increases the slab overload rapidly (Section 1). Here, the interpretation is made focusing on necessary conditions for natural avalanches, not direct estimation of avalanche frequency and magnitude (see also Section 6.3). Notably, a warmer climate would enhance the sintering process and would reduce the age of the weak layer (Section 1) commonly in Asahikawa, Tokamachi and Kiso (Fig. 12), which decreased the probability of weak layer formation over northern Japan (Fig. 6). Hence, for the weak layer formation, an increase in air temperature due to the +4°C warmer climate robustly overcame the effect of the snowfall intensity change. Moreover, the reduction of weak layer age reduced the total amount of slab overload in many areas (Fig. 9) following the correlation between the overload and weak layer age (Fig. 12). However, the overload would increase in central Hokkaido (or Asahikawa) for PP/DF and FC/DH weak layers and in Hokuriku (or Tokamachi) for the FC/DH weak layer (Figs 9, 11) because the increasing speed of slab overload would be faster (Figs 12a, c) following the increase in heavy snowfall (Fig. 2a). In central Honshu (or Kiso), while the heavy snowfall in NDJFM would decrease (Fig. 2b), the increasing speed of the overload would be faster (Fig. 12e). Since this faster increase in the overload was direct evidence of the increase in the snowfall amount during the weak layer existence, the increase in the overload for PP/DF weak layer in central Honshu (or Kiso) (Figs 9ac, 11e) would be due to the change in snowfall intensity. Hence, for the slab overloads, the effect of the change in snowfall intensity overcame the effect of the increase in air temperature, depending on the area.

The spatial dependency of the increase in slab overload as described above would be partially related to the dependency of the heavy snowfall on climatological air temperature (or elevation) as well as Le Roux and others (2021), who demonstrated the dependency of the historical frequency of heavy snowfall on elevations in France. According to Kawase and others (2020), who performed a numerical experiment using an RCM of 1 km spatial resolution forced with d4PDF in Hokuriku, the frequency of heavy snowfall would be increased/decreased in the elevations above/below 1000 m, where approximately corresponded to −7°C climatological air temperature in winter. This dividing temperature (or elevation) was in line with our result that the slab overload for FC/DH would be increased in the mountainous region of Hokuriku (Figs 1, 9, 10). However, the dividing temperature was slightly warmer than the general value in other countries (O'Gorman and others, 2014) probably because of the strong dependency of snowfall on SST in Hokuriku (Kawase and others, 2013, 2020).

From a similar perspective of the contradictory effect of changes in air temperature and snowfall intensity on the avalanche potential, we could interpret the spatial distribution of the weak layer formation (Fig. 6) and that of slab overload (Fig. 9) in the historical experiment. The SI value for both PP/DF and FC/DH was relatively low in Asahikawa (Fig. 8), where the cold air temperature delays the increase of the snowpack's strength via the sintering process (Section 2). Although this slow increase in the strength might make the amount of slab overload larger, the amount in Asahikawa would not be so much as that in Tokamachi (Fig. 11) because of the slower increasing speed of overloads following the less snowfall intensity in Asahikawa (Figs 12a, b). Moreover, the SI value for PP/DF was often 1.5 (Fig. 8a), a threshold value of the weak layer, in Tokamachi, where there are frequent snowfall events (Section 2; Fig. 2b), whereas the SI value was relatively large in Kiso, which has little snowfall (Section 2; Fig. 2c). Similarly, owing to the frequent snowfall intensity in Tokamachi, the slab overload for both PP/DF and FC/DH was larger than that in other regions (Fig. 11). For FC/DH, the amount of slab overload in Tokamachi would be mostly double that for PP/DF (Subsection 5.4) because of a combination of longer weak layer age and stronger intensity of snowfall.

Notably, the projected slab overload is based on a set of rare cases because the probability of weak layer formation was <1% of the total SCD (Fig. 6). For example, in Kiso, in the +4°C experiment, the weak layer of PP/DF would exist for 0.3 d per season, as determined by multiplying the probability of weak layer formation (1%; Fig. 6c) and total SCD (30 d; Fig. 4). Considering that the weak layer of PP/DF would last ~3 d in Kiso (Fig. 12e), a set of weak layers that last for 3 d would be formed once every 10 years, as calculated by dividing the lifetime of the weak layer by the number of days per season the weak layer exists (0.3 d a⁻¹). The slab overload corresponding to these rare weak layers would be projected as the most frequent value of ~50 kg m⁻² in Kiso in the +4°C experiment (Fig. 11e). Similarly, in Hokuriku (or Tokamachi), where the slab overload for FC/DH would increase by ~30% (Fig. 9d), a return period of a set of weak layers would be 10 years based on the probability of weak layers (1%; Fig. 6d), the number of SCD (60 d) and a lifetime of the weak layer (6 d; Fig. 12d). Thus, the slab overload would be large once the weak layer formed and lasted for several days, whereas the return period of weak layer formation would be generally long.

In this study, the weak layer and associated slab overload were simply identified by seeking a minimum SI of <1.5 across the snowpack layers (Subsection 4.3). In this simple strategy, although the weak layer would be only one of the necessary conditions for avalanches, many previous studies have reported that the SI negatively correlated with avalanche activities in many countries (Jamieson and Johnston, 1998; Zeidler, 2004; Schweizer and others, 2006; Castebrunet and others, 2012; Richter and others, 2020; Techel and others, 2020, Reuter and others, 2022), including Japan (Nishimura and others, 2005; Hirashima and others, 2008; Abe and Hirashima, 2015; Hirashima, 2019). The total number of avalanches evaluated was more than thousands including 33 avalanches in Japan. Hence, the result of a decrease in the probability of weak layer formation indicated the potential decrease of avalanche formation. Even though it is very difficult to convert the probability of weak layer formation to the actual probability of avalanche formation (Subsection 6.3), we may convert the probability to that of avalanche danger level provided by forecasters. According to Techel and others (2020), ~65% of cases where the poor/very poor snowpack stability measured by the rutschblock test are comparable to the avalanche danger level ‘considerable/high’ in Switzerland. Applying this proportion to our cases simply, 65% of the probability of weak layer formation (Fig. 6) might be the probability of avalanche danger level ‘considerable/high’. Moreover, because the weak layer is mechanically the most frangible, the snow layers above the weak layer are highly expected to be released by a failure. Considering that the magnitude of the avalanches should correlate with slab overloads and the area of the avalanche release zone, the result of an increase in the slab overload indicated the potential increase in avalanche magnitude. Moreover, for dry slab avalanches, because an initial failure could propagate far with denser (harder) slabs (Gaume and others, 2015), the denser slabs (Fig. 13) might indicate a larger release zone.

The potential decrease in avalanche formation over northern Japan (Fig. 6) also agrees with the decreasing trend of avalanche occurrence in the recent decades in the French Alps (Eckert and others, 2013), Swiss Alps (Teich and others, 2012a) and the US Northern Rocky Mountains (Peitzsch and others, 2021) and the estimation based on the combination of a climate model, a physical snowpack model and past avalanche records in the French Alps (Castebrunet and others, 2014). The warmer climate perhaps reduces the avalanche commonly in many countries, including Japan. Moreover, the potential magnitude of avalanches would be projected to decrease in many areas in this study (Fig. 9), which agreed with the recent temporal trend in the French Alps (Eckert and others, 2013). Meanwhile, in this study, we also realized a potential increase in avalanche magnitude associated with FC/DH weak layers in the mountainous areas of Hokuriku and central Hokkaido and that associated with PP/DF weak layers in central Honshu (Fig. 9). This partially contrasting result to Eckert and others (2013) is probably because their statistical process with a massif scale (~200 km) omitted the change in smaller areas. In fact, from the perspective of a bigger spatial scale, the future change in slab overload was blurred in this study as well (Fig. 10). Ballesteros-Cánovas and others (2018), who addressed the long-term avalanche record dendrogeomorphically retrieved from tree-rings in the Himalaya, also made a similar conclusion about the increasing avalanche frequency and magnitude, but their result was probably due to an increase in wet avalanches, which was not considered in this study.

6.3. Limitations and future tasks

In this study, we estimated the potential effect due to global warming including the consideration of changes in heavy snowfall. However, the heavy snowfall in the mountainous regions of Hokuriku and Hokkaido was overestimated by the NHRCM because of the insufficient resolution of topography (Subsection 3.2; Kawase and others, 2016). Because of this problem, the probability of weak layer formation estimated in this study could have a positive bias. In Tokamachi, for example, the slab overload for PP/DF was increased by 100 kg m⁻² per 5 d in the historical experiment (Fig. 12c), indicating that the daily snowfall amount was 20 kg m⁻². Because the frequency of this snowfall amount of the NHRCM was five times that of the observation around Tokamachi (Kawase and others, 2016), the probability of weak layer formation would have a 500% bias. However, the result of decreasing the probability of weak layer formation responding to global warming would be robust because the probability would decrease over northern Japan without exception (Fig. 6). Moreover, a significant increase in the frequency of heavy snowfall would be projected commonly by both 20 and 5 km RCMs (Kawase and others, 2016, 2020; Inatsu and others, 2021).

We estimated avalanche potential using the SNOWPACK model based on the assumption that each gridpoint represented flat terrain without vegetation, although avalanches might strongly depend on topography, snow redistribution and forests (Bebi and others, 2009; Bühler and others, 2013; CAA, 2016). In nature, local wind influences snow deposition in complex terrain (Liston and Strum, 1998). Incoming longwave and shortwave radiations depend on the terrain because of the interactions between snowpacks' surfaces and radiations (Lehning and others, 2006) and sometimes enhance FC/DH formation. These local effects probably reduce SI and increase slab overload in local scales <20 km simulation of this study. Hence, the estimation of avalanche potential based on meteorological data with a higher resolution than that used in this study considering the local effects on snowpacks could show a broader shape of the PDF of SI and slab overload, especially in mountainous regions. However, numerical calculations with higher resolution require a large computational cost, and so, they cannot be practically performed at the current level of computer performance. Fortunately, this study is expected to at least provide the avalanche potential in the 20 km spatial scale because snow particles rarely move beyond the 20 km grid spacing of the NHRCM and FC/DH would typically decrease because of global warming over northern Japan (Fig. 5). It is also mentionable that the changes between the historical and +4°C experiments would not be affected by the local effects under an assumption of constant topography and vegetation.

As another limitation of this study, the potential estimated using the SI just indicates a necessary condition of avalanche (Section 1). Therefore, it is ideally better to directly estimate the probability of avalanche formation and associated magnitude based on the full consideration of the mechanical process for avalanche formation: (i) the failure initiation of the snow layer, (ii) crack propagation through slab and (iii) slab tensile support (Reuter and Schweizer, 2018). Given this scientific trend, the critical crack length (CCL), the length of a crack propagating spontaneously after the initial failure of the snowpack's layer, is recently proposed as a new index to diagnose the mechanical process (ii) for a dry snow slab avalanche associated with the FC/DH weak layer (Gaume and Reuter, 2017; Gaume and others, 2017; Richter and others, 2019). Hence, the usage of CCL may relieve the problem that the SI just indicates a necessary condition. However, at present, it is difficult to define the weak layer using CCL properly due to lack of scientific work addressing the relationship between CCL and natural snowpack instability in Japan (e.g. Reuter and others, 2015; Gaume and Reuter, 2017; Rosendahl and Weißgraeber, 2020). As another approach bypassing the full consideration of the mechanical process, one may exploit the statistical relationship between avalanche and snowpack stability (Castebrunet and others, 2014). However, it requires enormous past avalanche datasets, which are not yet constructed in Japan, so it was also difficult to apply the procedure. Therefore, for further improvement of this study, both a huge dataset of past avalanches and a numerical model with the full consideration of the mechanical process are required.

6.4. Potential implications

The findings of this study indicated a decrease in avalanche occurrence and an increase in avalanche magnitude depending on areas, which potentially affect political/civilian adaptation to avalanche hazards. In several countries, including Japan, to mitigate avalanche hazard risk, people develop and maintain forest and artificial structures (Weir, 2002; Schneebeli and Bebi, 2004; Brang and others, 2006; Höller, 2007) to prevent avalanche formation and to decelerate avalanche flow (Kamiishi, 2004; Aiura, 2005; Höller, 2007; Bebi and others, 2009; Teich and others, 2012b; Takeuchi and others, 2018). These silvicultural/engineering managements are often implemented in areas where avalanches have occurred in the past or where high avalanche potential has been identified (Bebi and others, 2009; Bühler and others, 2013; CAA, 2016). Owing to the increase in avalanche magnitude, it might be necessary to reinforce the measures, which require economic costs. Meanwhile, fortunately, our results indicated a decrease in avalanche occurrence and magnitude in many areas, so the economic cost of avalanche management might be reduced in the future. Therefore, we expect to adapt to the changes in avalanches if the current measures were optimized through risk assessments improved by solving the limitations listed earlier (Subsection 6.3).

7. Conclusions

In this study, we predicted the effect of global warming on avalanche potential by forcing the SNOWPACK model with large-ensemble climate model results from d4PDF. The numerical snowpack calculation was performed for 1800 winter simulations in the historical climate (1951–2010) and for 1800 winter simulations in the future climate, of which the latter assumed that the global mean air temperature will increase by 4°C from the preindustrial condition. The SM-SWE reproduced in the historical experiment agreed with that of the long-term observation at Tokamachi from 1939 to 2020 (Fig. 3c). Although the model had problems of insufficient consideration of topography, vegetation and snow redistribution that potentially biased the calculation result, the results showed that global warming would decrease the SM-SWE over northern Japan (Fig. 3b), decrease SCD overall (Fig. 4b), increase PP/DF fraction (Figs 5a, c) and decrease FC/DH fraction (Figs 5b, d), which agreed well with previous studies. The avalanche potential was indicated by the weak layers and associated slab overloads, identified by the SI. Global warming would decrease the probability of weak layer formation over northern Japan independent of weak layer types, such as PP/DF and FC/DH (Fig. 6), indicating a potential decrease in avalanche formation. However, the slab overload would increase for the PP/DF weak layer in central Hokkaido and central Honshu and would increase for the FC/DH weak layer in the mountainous area of Hokuriku and central Hokkaido (Fig. 9), although the age of the weak layer would decrease by a few days over northern Japan (Fig. 12). This was because the effect of the increase in heavy snowfall, which increases the slab overload, overcame the effect of a shortened weak layer lifetime due to a temperature increase, which decreases the slab overload (Fig. 12). Because the slab above the weak layer could be released by a failure, an increase in slab overload indicated a potential increase in avalanche magnitude.