Journal of Glaciology

Large-ensemble climate simulations to assesschanges in snow stability over northern Japan
Yuta Katsuyama

ArticleCite this article: Katsuyama Y, Katsushima T,
Takeuchi Y (2022). Large-ensemble climate
simulations to assess changes in snow stability
over northern Japan. Journal of Glaciology
1–14. https://doi.org/10.1017/jog.2022.85
Received: 12 July 2021
Revised: 5 September 2022
Accepted: 9 September 2022
Key words:
Avalanches; climate change; snow
metamorphosis
Author for correspondence:
Yuta Katsuyama,
E-mail: ykatsuyama2020@affrc.go.jp, Takafumi Katsushima and Yukari Takeuchi

Tohkamachi Experimental Station, Forestry and Forest Products Research Institute, Tokamachi 9480013, Japan

Abstract
To examine the influence of global warming including increased heavy snowfall frequency on thepotential 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 thenatural stability index, would decrease all over the area because of the shorter age of the weaklayers caused by a warmer climate, indirectly indicating a potential decrease in avalanche frequency. However, because of increased heavy snowfall frequency, slab overload would increaseby 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© The Author(s), 2022. Published byCambridge University Press on behalf of The
International Glaciological Society. This is anOpen Access article, distributed under theterms of the Creative Commons Attribution
licence (http://creativecommons.org/licenses/
by/4.0/), which permits unrestricted re-use,
distribution and reproduction, provided theoriginal article is properly cited.

cambridge.org/jog

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 andothers, 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 inrecent decades below an elevation of 1200 m in the Vosges Mountains in France (Giaconaand others, 2021). Nevertheless, an increase in wet avalanches would be expected becauseof 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 andothers, 2021). Hence, the spatial variation of avalanche formation and magnitude might beinfluenced 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) calculatedusing a physical snowpack model imposed with perturbed atmospheric data. Then, theyfound 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 theincrease 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 meteorologicalstate of the present climate and by seeking a similar meteorological state between the present and future climates. Then, they determined that avalanche activities decreased andhttps://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

2

increased at low and high altitudes, respectively. However, theirstatistical 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 atpresent.
Notably, in the assessment, increases in air temperature andheavy 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 weaklayers. This effect of temperature increase would also reduce theamount of overlying slab above the weak layer because the slab
was fed by less amount of snowfall during the shorter lifetimeof 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 snowpackstability, researchers have employed physical-based snowpack
models, such as SNOWPACK (Bartelt and Lehning, 2002) andCrocus (Brun and others, 1992), which produce a time evolutionof 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 forpolicy decision making for future climate change (d4PDF). Thisdataset 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 increasein mountainous areas in central Japan and would decrease incoastal 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 thatused by Katsuyama and others (2020), who estimated the effectof global warming on the amount and grain types of snowpackover 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 andshear 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 thehttps://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

Yuta Katsuyama and others

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 theirextremes given by d4PDF. Then, the snow water equivalent(SWE), snow-covered days (SCD) and snow grain types calculatedusing 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 featuressouthwest–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 thevapor produced by the Sea of Japan (e.g. Takano and others,2008). In particular, in the Hokuriku region (Fig. 1), snowfall isheavier than that in other regions because of an updraft enhancedby the mountains (Houze, 2014) and moisture flux convergencein 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.
Asahikawa, Tokamachi and Kiso were selected as representative locations of central Hokkaido, Hokuriku and central
Honshu, respectively (Fig. 1b). Asahikawa is characterized ashaving 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 mainlycomprises dry snow and sometimes faceted crystals and depthhoar (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,

Journal of Glaciology

a

3

b

Fig. 1. (a) Topography of the Meteorological ResearchInstitute 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.

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 theprobability 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 usedinstead of data from manual observation. The PDF of SM-SWE
was constructed from a series of observed SM-SWEs.
3.2. d4PDF
In this study, the SNOWPACK model was forced with the resultsof the RCM in d4PDF (Mizuta and others, 2017). These data werecreated by dynamically downscaling using a nonhydrostatic RCM(NHRCM) with a 20 km spatial resolution from version 3.2 ofthe 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. theincrease in the global mean air temperature from the historicallevel was +3.55°C). In this +4°C experiment, 15 ensemble simulations were performed using different initial conditions andperturbated 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 fiveensemble 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 theperiod from 2051 to 2111 (Mizuta and others, 2017). Moreover,it should be noted that, in the ensemble experiment like thisstudy, the hourly (or daily) time series are not comparable tothe observations, but are comparable in terms of distribution(or associated statistics such as mean and std dev.).
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 andothers, 2016), including snowfall that occurs when the air temperature is ∼0°C, which is a characteristic of Hokuriku (Ohbaand 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 andoverestimated 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
Table 2. Selected d4PDF
Historical experiment
PeriodTable 1. Observation intervals of SWE for each period at Tokamachi, Niigata,Japan

EnsemblememberSea surface
temperature

PeriodTotal

1939/40–1946/47
1947/48
1948/49–1951/52
1952/53–1956/57
1957/58
1958/59–2020/21

https://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

Interval [days]
5
10
5
1
10
1

a

1 August 1951 to 31
July 2011
m001–m010, m021–
m030, m041–m050
Observation

1800 seasons+4°C experiment1 August 2051 to 31 July 2111
m101–m105
Projection from CCSM4a, GFDL
CM3b, HadGEM2c, MIROC5d,
MPI-ESMe and MRI-CGCM3f
1800 seasons

Community Climate System Model version 4 (CCSM4).
Geophysical Fluid Dynamics Laboratory Climate Model version 3(GFDL CM3).
Hadley Centre Global Environment Model version 2 (HadGEM2).
dModel 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 Modelversion 3 (MPI-CGCM3).
bc

4provide complete hourly time series of snowpacks’ structures andphysical properties, such as grain type, density and shear strength(Bartelt and Lehning, 2002), for all ensemble members of both thehistorical and +4°C experiments (Supplement 1 provides someexample time series of SWE at Tokamachi). The calculation
period for the historical experiment was from 1951/1952 winterto 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 andstopped on July 31 of the following year. The SNOWPACK
model was forced with seven atmospheric variables (air temperature, wind speed, relative humidity, incoming shortwave andlongwave radiations, and solid and liquid precipitations) every
hour. The model was initialized with 10 m of soil layers with a + 2.3 [◦ C]°C, where T
 is the climatology
temperature of 0.97Tof 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 distributionof snowpacks simply by looping the model at each gridpoint. The
grid configuration and elevation at each gridpoint were the sameas 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 thisassumption will be discussed in Section 6. See Katsuyama andothers (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 thenumerical 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 andothers, 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 inJapan, the PDF of SM-SWE of large-ensemble simulations wascompared 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 constructedfrom 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 asmean and std dev.). Therefore, for example, the daily time seriesof snow depth in a specific year were not correlated between thecalculation 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 wasmore 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

https://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

Yuta Katsuyama and others

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 andwinter is a weighted mean, x = i fi di / i di , 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 historicaland +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 anindex 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 toestimate the potential of avalanches, assuming a uniform slope(Föhn, 1987; Jamieson and Johnston, 1998):
Su =

s
,
w cos u sin u(1)where S is the SI, w [Pa] is the vertical stress due to the weight ofthe 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 mechanicallyweak 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, S38, 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 thestructural 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 weaklayer 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 bykinetic 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 theminimum 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 unstablesnowpacks (Perla, 1977; Jamieson and Johnston, 1995; Lehningand others, 2004; Hirashima and others, 2008). This thresholdmethod using SI estimates the best probability of capturing anecessary 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 weaklayers of PP/DF and FC/DH are independent. Finally, the potential magnitude of avalanches was identified on the basis of theamount of the slab overlying the weak layer. The weak layer

Journal of Glaciology

5and 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−3, 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:

P|S38 ,1.5 =

Nm

i=1

24Ns ij=1

H0 (1.5 − S38i,j )
24Nsi

,

a

b

(2)

where Nm (=1800) is the total number of winters of each experiment, Ns is the number of SCD and H0 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 forall ensemble members in each of the historical and +4°C experiments. In the analysis for the three locations, we selected the threenearest grid points on the model to a location. Then, the PDFs of
SI, slab overload and slab density at the locations were constructedusing 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 fromthe 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).

cFig. 2. Cumulative frequency distribution of daily snowfall amount from November toMarch in the following year at (a) Asahikawa, (b) Tokamachi and (c) Kiso. The blackand 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
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 thehistorical 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−2 d−1 (Fig. 2a). This pattern ofchange in CFD directly means the increase in the frequency ofheavy snowfall. However, even with the increases in the heavy
snowfall, the mean daily precipitation in NDJFM would notchange (figure not shown). In Tokamachi, the snowfall amountwith a longer/shorter return period than once per 5 winters,that was ∼75 kg m−2 d−1, 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% comparingthat 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 increasein 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).

https://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

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 areaof Hokuriku because of the heavy snowfall (Section 2; Fig. 2b).
In Hokkaido, the SWE was ∼600 kg m−2 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−2 inthe 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 thehistorical experiment. In mountainous areas, the decrease would
be relatively small, i.e. ∼20–30%.
The PDF of the SM-SWE was also compared with that ofTokamachi (Fig. 3c). It showed a broad distribution with amean of 800 kg m−2 and a range of 100–1500 kg m−2. 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−2,
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−2) was less than the mean value in the historicalexperiment. Notably, the double-peak shape of the PDF in thehistorical experiment was simply because multiple grid points
on the model were combined (Subsection 4.1).

6

Yuta Katsuyama and others

a

b

c

Fig. 3. (a) Mean of the seasonal maximum snow water equivalent (SM-SWE) [kg m−2] 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.

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, andthe 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, whichmostly indicated no snow-covered.
Figure 5 shows the thickness fraction of grain types. The fraction of PP/DF largely increased from ∼20% in the historicalexperiment 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 toexist 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).

a

5.3. Weak layer formation
Figure 6 shows the existence probabilities of weak layers where theminimum 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

a

b

c

d

b

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

https://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

Fig. 5. (a) Thickness fraction of precipitation particles and decomposing and fragmented precipitation particles (PP/DF) and (b) that of faceted crystals and depthhoar (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−2 in the historical experiment, were masked out (Fig. 3a).

Journal of Glaciology

7

a
a

b

b

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

cdFig. 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/DFand (d) FC/DH in the +4°C experiment. Areas, where the SM-SWE was <50 kg m−2 in thehistorical experiment, were masked out (Fig. 3a).

corresponding to the spatial distribution of SWE (Fig. 3a): theprobability of weak layer formation was more than 20% in themountainous region of Hokuriku and ∼10% in Hokkaido.
Following this spatial pattern, the probability was well correlatedwith 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).
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 anormal 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, aathreshold 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/DFand 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/DFor FC/DH is not always included in the snowpack for every
time step.
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 themountainous area of Hokuriku and central Hokkaido in the+4°C experiment (Fig. 6c). For FC/DH, in all regions, the weaklayer 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 distributionof 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 historicalexperiment, the mean overload above the weak layer of PP/DF
was ∼100 kg m−2 in the mountainous region of Hokuriku,
which was considerably larger than that of the other regions,

a

b

c

d

b

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

https://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

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 theslab overload of the historical experiment. Areas, where the SM-SWE was <50 kg m−2
in the historical experiment, were masked out (Fig. 3a).

8

Yuta Katsuyama and othersabFig. 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 theseasonal maximum SWE was <50 kg m−2 in the historical experiment, were excluded(Fig. 3a). The bins for elevation are 100 m intervals.

∼40 kg m−2 in Hokkaido and the central Honshu (Fig. 9a). The
spatial distribution of overload for FC/DH was similar to thatfor 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 ofoverload 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, theamount of overload for FC/DH would be ∼160 kg m−2, much larger than that for PP/DF. However, in Hokkaido and central
Honshu, the overload for FC/DH in the historical experimentwas ∼40 kg m−2, mostly the same amount as that for PP/DF(Figs 9a, b). Following this spatial pattern, the overloads forboth PP/DF and FC/DH were well correlated with elevation(Fig. 10).
Figure 11 shows the PDF of the slab overload in the three locations. The overload for PP/DF in the historical experiment was

a

b

c

d

e

f

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 theresults of the historical and +4°C experiments, respectively.

https://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

the most frequent at ∼50 kg m−2 in Asahikawa (Fig. 11a) andTokamachi (Fig. 11c) and at ∼20 kg m−2 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−2 inAsahikawa, Tokamachi and Kiso, respectively. Similarly, theprojected 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−2, the most frequent value being 150 kg m−2(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 increasein the mountainous region of Hokuriku and central Hokkaido for
the FC/DH weak layer (Fig. 9d), whereas that would decrease inmany other areas (Fig. 9). In central Honshu, in particular, themean slab overload for PP/DF would increase by ∼30% relative
to the historical experiment (Fig. 9c). Because the PDF peak ofslab overload for PP/DF would significantly shift from 10 to
50 kg m−2 in Kiso (Fig. 11e), the increase in the mean slab
overload in central Honshu would occur simply because of thepeak shift. In Asahikawa, the probability of slab overload for
20 kg m−2 in the historical experiment would be half in the+4°C experiments (Fig. 11a), and this would cause a slightincrease 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 aslight increase would be projected in Hokuriku and central
Hokkaido. This increase in mean slab overload would be due toan increase in probability corresponding to a 180 kg m−2 slab
overload in Tokamachi (Fig. 11d) and that corresponding to a
100 kg m−2 slab overload in Asahikawa (Fig. 11b) in the PDFs
between the historical and +4°C experiments. In terms of thedependency 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 andslab 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 historicalexperiment. 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 threelocations, the slab overloads increased with an increase in weaklayer 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 thetwo sides of the same coin via the settlement process. In
Asahikawa, following the increase in the probability of the overload for 50 kg m−2 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−3 compared to the historicalexperiment (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−3, the probability which the avalanche release zone

Journal of Glaciology

9

a

b

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

c

d

e

f

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−3 and 10−4, respectively.

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−2 for the PP/DF (Fig. 11e), the PDF of slab density

a

b

c

d

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 withthe previous results. The SWE spatial pattern in the historicalexperiment (Fig. 3a) was similar to that found by Hosaka andothers (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−2 in the coastal region along the Sea ofJapan 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−2 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 andothers, 2005; Katsuyama and others, 2020) because the amountof snowpack at a high altitude was correlated with precipitation
rather than air temperature (Sun and others, 2016). The samedecreasing pattern was obtained in this study. The monotonicaldecrease in SCD by ∼30 d in Hokkaido agreed well with the resultof 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 estimationof primary grain type over Japan based on a statistical relationshipbetween atmospheric climatology and grain types observed (Inoueand Yokoyama, 2003; Ishizaka, 2008).
6.2. Weak layer formation and associated slab overloads anddensityefFig. 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 thehistorical and +4°C experiment, respectively. The black and red text beside thegray vertical dashed line shows the probability integrated over the density larger
than 150 kg m−3 in the historical and +4°C experiment, respectively.

https://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

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 theincrease in precipitation intensity weakens the snowpack andincreases the slab overload rapidly (Section 1). Here, the interpretation is made focusing on necessary conditions for natural
avalanches, not direct estimation of avalanche frequency andmagnitude (see also Section 6.3). Notably, a warmer climate
would enhance the sintering process and would reduce the ageof the weak layer (Section 1) commonly in Asahikawa,
Tokamachi and Kiso (Fig. 12), which decreased the probabilityof weak layer formation over northern Japan (Fig. 6). Hence,for the weak layer formation, an increase in air temperature dueto the +4°C warmer climate robustly overcame the effect of the10snowfall 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 andweak layer age (Fig. 12). However, the overload would increasein 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 fasterincrease in the overload was direct evidence of the increase in thesnowfall 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 ofthe heavy snowfall on climatological air temperature (or elevation)
as well as Le Roux and others (2021), who demonstrated thedependency 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 ofheavy 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 regionof 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 ofchanges 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 thisslow increase in the strength might make the amount of slab overload larger, the amount in Asahikawa would not be so much asthat in Tokamachi (Fig. 11) because of the slower increasing
speed of overloads following the less snowfall intensity inAsahikawa (Figs 12a, b). Moreover, the SI value for PP/DF wasoften 1.5 (Fig. 8a), a threshold value of the weak layer, inTokamachi, 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 ofsnowfall.
Notably, the projected slab overload is based on a set of rarecases 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 thatthe weak layer of PP/DF would last ∼3 d in Kiso (Fig. 12e), a

https://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

Yuta Katsuyama and others

set of weak layers that last for 3 d would be formed once every
10 years, as calculated by dividing the lifetime of the weaklayer by the number of days per season the weak layer exists(0.3 d a−1). The slab overload corresponding to these rare
weak layers would be projected as the most frequent value of∼50 kg m−2 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 periodof a set of weak layers would be 10 years based on the probabilityof 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 thesnowpack layers (Subsection 4.3). In this simple strategy, although
the weak layer would be only one of the necessary conditions foravalanches, many previous studies have reported that the SI negatively correlated with avalanche activities in many countries(Jamieson and Johnston, 1998; Zeidler, 2004; Schweizer andothers, 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 theprobability 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 byforecasters. According to Techel and others (2020), ∼65% of cases
where the poor/very poor snowpack stability measured by therutschblock 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 themagnitude of the avalanches should correlate with slab overloads
and the area of the avalanche release zone, the result of anincrease in the slab overload indicated the potential increase inavalanche 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 andothers, 2013), Swiss Alps (Teich and others, 2012a) and the US
Northern Rocky Mountains (Peitzsch and others, 2021) and theestimation 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 andcentral 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

Journal of Glaciology

process with a massif scale (∼200 km) omitted the change insmaller areas. In fact, from the perspective of a bigger spatial
scale, the future change in slab overload was blurred in thisstudy 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 asimilar 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 havea positive bias. In Tokamachi, for example, the slab overload for
PP/DF was increased by 100 kg m−2 per 5 d in the historicalexperiment (Fig. 12c), indicating that the daily snowfall amountwas 20 kg m−2. Because the frequency of this snowfall amountof the NHRCM was five times that of the observation around
Tokamachi (Kawase and others, 2016), the probability of weaklayer 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 mightstrongly 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 ofavalanche 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 andslab 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, thisstudy 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 andvegetation.
As another limitation of this study, the potential estimatedusing the SI just indicates a necessary condition of avalanche(Section 1). Therefore, it is ideally better to directly estimate theprobability 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

https://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

11(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, isrecently 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 inJapan (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 constructedin Japan, so it was also difficult to apply the procedure. Therefore,
for further improvement of this study, both a huge dataset of pastavalanches 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 toavalanche 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 adecrease in avalanche occurrence and magnitude in many areas,so the economic cost of avalanche management might be reducedin the future. Therefore, we expect to adapt to the changes inavalanches 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 withlarge-ensemble climate model results from d4PDF. The numericalsnowpack calculation was performed for 1800 winter simulationsin 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 historicalexperiment 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) anddecrease FC/DH fraction (Figs 5b, d), which agreed well with previous studies. The avalanche potential was indicated by the weak

12layers 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 themountainous 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 effectof 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.
Supplementary material. The supplementary material for this article canbe found at https://doi.org/10.1017/jog.2022.85.
Acknowledgements. We thank an editor and three anonymous reviewers
for giving us insightful comments that helped us improve the original manuscript. This study was supported by the JSPS KAKENHI grant (No. 20K22437
and 22K14458) and research grant of the Forestry and Forest Products
Research Institute (No. 202003). We used the supercomputer of AFFRIT,
MAFF, Japan and d4PDF produced with the Earth Simulator jointly by science
programs (SOUSEI, TOUGOU, SI-CAT and DIAS) of the Ministry of
Education, Culture, Sports, Science and Technology (MEXT), Japan. A python
library of mpi4py was partly used to calculate some statistics of numericalsimulation results (Dalcín and others, 2008).

References
Abe O and Hirashima H (2015) Evaluation of an avalanche forecasting system
based on the development process of faceted crystals. Seppyo 77(1), 37–45,(in Japanese).
Aiura H (2005) Stand density to control slope snow cover movements. Journalof the Japanese Society 87(1), 73–79, (in Japanese).
Ballesteros-Cánovas JA, Trappmann D, Madrigal-González J, Eckert N and
Stoffel M (2018) Climate warming enhances snow avalanche risk in theWestern Himalayas. Proceedings of the National Academy of Sciences of
the USA 115(13), 3410–3415. doi: 10.1073/pnas.1716913115
Bartelt P and Lehning M (2002) A physical SNOWPACK model for the Swissavalanche warning Part I: numerical model. Cold Regions Science andTechnology 35(3), 123–145. doi: 10.1016/S0165-232X(02)00074-5
Bebi P, Kulakowski D and Rixen C (2009) Snow avalanche disturbances inforest ecosystems-state of research and implications for management.
Forest Ecology and Management 257(9), 1883–1892. doi: 10.1016/j.foreco.
2009.01.050
Birkeland KW, Johnson RF and Schmidt DS (1998) Near-surface faceted
crystals formed by diurnal recrystallization: a case study of weak layer formation in the mountain snowpack and its contribution to snow avalanches.
Arctic, Antarctic, and Alpine Research 30(2), 200–204. doi: doi:10.2307/
1552135
Brang P and 5 others (2006) Management of protection forests in theEuropean Alps: an overview. Forest Snow and Landscape Research 80(1),23–44.
Brun E, David P, Sudul M and Brunot G (1992) A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting. Journalof Glaciology 38(128), 13–22.
Bühler Y and 5 others 2013) Automated identification of potential snow avalanche release areas based on digital elevation models. Natural Hazards andEarth System Sciences 13(5), 1321–1335. doi: 10.5194/nhess-13-1321-2013
CAA (2016) Technical Aspects of Snow Avalanche Risk Management─Resourcesand Guidelines for Avalanche Practitioners in Canada. Revelstoke, British
Columbia: Canadian Avalanche Association (CAA).
Castebrunet H, Eckert N and Giraud G (2012) Snow and weather climatic
control on snow avalanche occurrence fluctuations over 50 yr in the
French Alps. Climate of the Past 8, 855–875. doi: 10.5194/cp-8-855-2012https://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

Yuta Katsuyama and others
Castebrunet H, Eckert N, Giraud G, Durand Y and Morin S (2014)
Projected changes of snow conditions and avalanche activity in a warming
climate: the French Alps over the 2020–2050 and 2070–2100 periods. The
Cryosphere 8, 1673–1697. doi: 10.5194/tc-8-1673-2014
Dalcín L, Paz R, Storti M and D’Elía J (2008) MPI for python: performance
improvements and MPI-2 extensions. The Journal of Parallel andDistributed Computing 68(5), 655–662. doi: https://doi.org/10.1016/j.jpdc.
2007.09.005.
Eckert N, Keylock CJ, Castebrunet H, Lavigne A and Naaim M (2013)
Temporal trends in avalanche activity in the French Alps and subregions:from occurrences and runout altitudes to unsteady return periods.
Journal of Glaciology 59(213), 93–114. doi: 10.3189/2013JoG12J091
Endo Y (1993) Forecasting of direct-action avalanches in terms of snow accumulation rates. Seppyo 55(2), 113–120. doi: 10.5331/seppyo.55.113 (inJapanese).
Föhn PMB (1987) The stability index and various triggering mechanisms.
International Association of Scientific Hydrology 162, 195–214.
Gaume J, Chambon G, Eckert N, Naaim M and Schweizer J (2015) Influenceof weak layer heterogeneity and slab properties on slab tensile failure propensity and avalanche release area. The Cryosphere 9, 795–804. doi: 10.
5194/tc-9-795-2015
Gaume J and Reuter B (2017) Assessing snow instability in skier-triggeredsnow slab avalanches by combining failure initiation and crack propagation.
Cold Regions Science and Technology 144, 6–15. doi: 10.1016/j.coldregions.
2017.05.011
Gaume J, van Herwijnen A, Chambon G, Wever N and Schweizer J (2017)
Snow fracture in relation to slab avalanche release: critical state for the onsetof crack propagation. The Cryosphere 11, 217–228. doi: 10.5194/tc-11-217-2017
Giacona F and 7 others (2021) Upslope migration of snow avalanches in a
warming climate. Proceedings of the National Academy of Sciences of the
USA 118(44), e2107306118. doi: 10.1073/pnas.2107306118
Hirashima H (2019) Numerical snowpack model simulation schemes for avalanche prediction in Japan. Bulletin of Glaciological Research 37S, 37–41.
doi: 10.5331/bgr.18SW02
Hirashima H, Abe O, Sato A and Lehning M (2009) An adjustment for
kinetic growth metamorphism to improve shear strength parameterizationin the SNOWPACK model. Cold Regions Science and Technology 59(2–3),
169–177. doi: 10.1016/j.coldregions.2009.05.001
Hirashima H, Nishimura K, Baba E, Hachikubo A and Lehning M (2004)
SNOWPACK model simulations for snow in Hokkaido, Japan. Annals of
Glaciology 38, 123–129.
Hirashima H, Nishimura K, Yamaguchi S, Sato A and Lehning M (2008)
Avalanche forecasting in a heavy snowfall area using the snowpack
model. Cold Regions Science and Technology 51(2–3), 191–203. doi: 10.
1016/j.coldregions.2007.05.013
Hirota T, Fukumoto M, Shirooka R and Muramatsu K (1995) Simple
method of estimating daily mean soil temperature by using the force-restore
model. Journal of Agricultural Meteorology 51(3), 269–277. doi: 10.2480/
agrmet.51.269, (in Japanese).
Hock R and 14 others (2019) High Mountain Areas. IPCC Special Report on
the Ocean and Cryosphere in a Changing Climate.
Höller P (2007) Avalanche hazards and mitigation in Austria: a review.
Natural Hazards 43(1), 81–101. doi: 10.1007/s11069-007-9109-2
Hosaka M, Nohara D and Kitoh A (2005) Changes in snow cover and snow
water equivalent due to global warming simulated by a 20 km-mesh global
atmospheric model. SOLA 1, 93–96. doi: 10.2151/sola.2005-025
Houze AR (2014) Chapter 12 - Clouds and precipitation associated with hills
and mountains. In Houze RA (ed.), Cloud Dynamics. San Diego, CA:Elsevier Academic Press. International Geophysics 104, 369–402. doi:
10.1016/B978-0-12-374266-7.00012-3
Inatsu M, Kawazoe S and Mori M (2021) Trends and prediction of heavy
snowfall in Hokkaido, Japan as an application of self-organizing map.
Journal of Applied Meteorology and Climatology 60(10), 1483–1494. doi:
10.1175/JAMC-D-21-0085.1
Inoue S and Yokoyama K (2003) Estimates of snowfall depth, maximum
snow depth, and snow pack environments under global warming inJapan from five sets of predicted data. Journal of Agricultural Meteorology
59(3), 227–236. doi: 10.2480/agrmet.59.227
Ishizaka M (2008) Reassessment of climatic conditions in ‘depth-hoar region’and new map for climatic division of snow-covered areas in Japan based on
the new conditions. Seppyo 70(1), 3–13. doi: doi.org/10.5331/seppyo.70.1_3,(in Japanese).

Journal of Glaciology
Jamieson JB and Johnston CD (1995) Shear frame stability parameters forlarge-scale avalanche forecasting. Annals of Glaciology 18(1), 268–273.
doi: 10.3189/S0260305500011630
Jamieson JB and Johnston CD (1998) Refinements to the stability index forskier-triggered dry-slab avalanches. Annals of Glaciology 26, 296–302. doi:
10.3189/1998AoG26-1-296-302
JMA (2021) Web page of Japan Meteorological Agency (JMA). Available athttps://www.jma.go.jp/jma/index.html.
Kamiishi I (2004) Counter-measure for road avalanche. Journal of SnowEngineering of Japan 20(2), 69–75, (in Japanese).
Katsuyama Y, Inatsu M, Nakamura K and Matoba S (2017) Global warming
response of snowpack at mountain range in northern Japan estimated using
multiple dynamically downscaled data. Cold Regions Science and Technology
136, 62–71. doi: 10.1016/j.coldregions.2017.01.006
Katsuyama Y, Inatsu M and Shirakawa T (2020) Response of snowpack to +2 °C
global warming in Hokkaido, Japan. Journal of Glaciology 66(255), 83–96.
doi: 10.1017/jog.2019.85
Kawase H and 6 others (2013) Altitude dependency of future snow coverchanges over Central Japan evaluated by as regional climate model.
Journal of Geophysical Research-Atmospheres 118, 12444–12457. doi: 10.
1002/2013JD020429
Kawase H and 6 others (2016) Enhancement of heavy daily snowfall in central
Japan due to global warming as projected by large ensemble of regional climate simulations. Climatic Change 139(2), 265–278. doi: 10.1007/
s10584-016-1781-3
Kawase H and 10 others (2020) Changes in extremely heavy and light snowcover winters due to global warming over high mountainous areas in central
Japan. Progress in Earth and Planetary Sciences 7, 10. doi: https://doi.org/10.
1186/s40645-020-0322-x.
Lavigne A, Eckert N, Bel L and Parent E (2015) Adding expert contributions
to the spatiotemporal modelling of avalanche activity under different climatic influences. Journal of the Royal Statistical Society C 64, 651–671.
doi: 10.1111/rssc.12095
Lehning M and 5 others (2006) ALPINE3D: a detailed model of mountainsurface processes and its application to snow hydrology. Hydrological
Processes 20(10), 2111–2128. doi: 10.1002/hyp.6204
Lehning M, Bartelt P, Brown B and Fierz C (2002) A physical SNOWPACK
model for the Swiss avalanche warning Part III : meteorological forcing,thin layer formation and evaluation. Cold Regions Science and Technology
35, 169–184. doi: 10.1016/s0165-232x(02)00074-5
Lehning M, Fierz C, Brown B and Jamieson B (2004) Modeling snow
instability with the snow-cover model SNOWPACK. Annals of Glaciology
38, 331–338. doi: 10.3189/172756404781815220
Le Roux E, Evin G, Eckert N, Blanchet J and Morin S (2021)
Elevation-dependent trends in extreme snowfall in the French Alps from1959 to 2019. The Cryosphere 15, 4335–4356. doi: 10.5194/tc-15-4335-2021
Liston EG and Strum M (1998) A snow-transport model for complex terrain.
Journal of Glaciology 44(148), 498–516.
López-Moreno JI, Pomeroy JW, Revuelto J and Vicente-Serrano SM (2012)
Response of snow processes to climate change: spatial variability in a small
basin in the Spanish Pyrenees. Hydrological Processes 27, 2637–2650. doi:
10.1002/hyp.9408
Martin E, Giraud G, Lejeune Y and Boudart G (2001) Impact of a climatechange on avalanche hazard. Annals of Glaciology 32, 163–167.
McClung D and Schaerer P (2006) The Avalanche Handbook, 3rd Edn.
Seattle, WA: Mountaineers Books.
Mitterer C, Hirashima H and Schweizer J (2011) Wet-snow instabilities:
comparison of measured and modelled liquid water content and snow
stratigraphy. Annals of Glaciology 52(58), 201–208. doi: 10.3189/
172756411797252077
Mizuta R and 30 others (2017) Over 5,000 years of ensemble future climate
simulations by 60-km global and 20-km regional atmospheric models.
Bulletin of the American Meteorological Society 98(7), 1383–1398. doi: 10.
1175/BAMS-D-16-0099.1
Monti F, Gaume J, Van Herwijnen A and Schweizer J (2016) Snow instability evaluation: calculating the skier-induced stress in a multi-layered snowpack. Natural Hazards and Earth System Sciences 16(3), 775–788. doi:
doi:10.5194/nhess-16-775-2016
Morin S and 16 others (2020) Application of physical snowpack models insupport of operational avalanche hazard forecasting: a status report on current implementations and prospects for the future. Cold Regions Science andTechnology 170, 102910. doi: 10.1016/j.coldregions.2019.102910

https://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

13
Nagata M (1987) On the structure of a convergent cloud band over the Japansea in winter; A prediction experiment. Journal of the Meteorological Societyof Japan 65(6), 871–883. doi: 10.2151/jmsj1965.65.6_871
Nishimura K, Baba E, Hirashima H and Lehning M (2005) Application of thesnow cover model SNOWPACK to snow avalanche warning in Niseko, Japan.
Cold Regions Science and Technology 43(1–2), 62–70. doi: 10.1016/j.coldregions.2005.05.007
O’Gorman P (2014) Contrasting responses of mean and extreme snowfall toclimate change. Nature 512, 416–418. doi: 10.1038/nature13625
Ohba M and Sugimoto S (2020) Impacts of climate change on heavy wet
snowfall in Japan. Climate Dynamics 54(5–6), 3151–3164. doi: 10.1007/
s00382-020-05163-z
Peitzsch EH, Pederson GT, Birkeland KW, Hendrikx J and Fagre DB (2021)
Climate drivers of large magnitude snow avalanche years in the U.S. northern Rocky Mountains. Scientific Reports 11, 10032. doi: 10.1038/
s41598-021-89547-z
Perla R (1977) Slab avalanche measurements. Canadian Geotechnical Journal
14(2), 206–213. doi: 10.1139/t77-021
Reuter B and 6 others (2022) Characterizing snow instability with avalanche
problem types derived from snow cover simulations. Cold Regions Scienceand Technology 194, 103462. doi: 10.1016/j.coldregions.2021.103462
Reuter B and Schweizer J (2018) Describing snow instability by failure initiation, crack propagation, and slab tensile support. Geophysical Research
Letters 45, 7019–7027. doi: 10.1029/2018GL078069
Reuter B, Schweizer J and van Herwijnen A (2015) A process-based
approach to estimate point snow instability. The Cryosphere 9, 837–847.
doi: 10.5194/tc-9-837-2015
Richter B, Schweizer J, Rotach MW and van Herwijnen A (2019) Validating
modeled critical crack length for crack propagation in the snow covermodel SNOWPACK. The Cryosphere 13, 3353–3366. doi: 10.5194/
tc-13-3353-2019
Richter B, van Herwijnen A, Rotach M W and Schweizer J (2020) Sensitivityof modeled snow stability data to meteorological input uncertainty. Natural
Hazards and Earth System Sciences 20, 2873–2888. doi: 10.5194/
nhess-20-2873-2020
Rosendahl PL and Weißgraeber P (2020) Modeling snow slab avalanchescaused by weak-layer failure – Part 2: coupled mixed-mode criterion forskier-triggered anticracks. The Cryosphere 14, 131–145.
Sasai T and 8 others (2019) Future projection of extreme heavy snowfall
events with a 5-km large ensemble regional climate simulation. Journal ofGeophysical Research-Atmospheres 124(24), 13975–13990. doi: 10.1029/
2019JD030781
Schneebeli M and Bebi P (2004) Snow and avalanche control. In Burley J, Evans J
and Youngquist J (eds), Encyclopedia of Forest Sciences. Amsterdam: Elsevier,pp. 397–402.
Schweizer J, Bellaire S, Fierz C, Lehning M and Pielmeier C (2006)
Evaluating and improving the stability predictions of the snow covermodel SNOWPACK. Cold Regions Science and Technology 46(1), 52–59.
doi: 10.1016/j.coldregions.2006.05.007
Strapazzon G and 5 others (2021) Effects of climate change on avalanche
accidents and survival. Frontiers in Physiology 12, 1–10. doi: 10.3389/
fphys.2021.639433
Sun F, Hall A, Schwartz M, Walton DB and Berg N (2016) Twenty-first-century
snowfall and snowpack changes over the Southern California mountains.
Journal of Climate 29(1), 91–110. doi: 10.1175/JCLI-D-15-0199.1
Takano Y, Tachibana Y and Iwamoto K (2008) Influences of large-scaleatmospheric circulation and local sea surface temperature on convective
activity over the sea of Japan in December. SOLA 4, 113–116. doi: 10.
2151/sola.2008-029
Takeuchi Y and 6 others (2019a) Meteorological statistics during the past 100
years from 1918 to 2017 at Tohkamachi Experimental Station, Forestry andForest Products Research Institute, Japan. Bulletin of Forestry and Forest
Products Research Institute 18(1), 35–99. doi: 10.20756/ffpri.18.1_35 (inJapanese).
Takeuchi Y, Katsushima T and Endo Y (2019b) Data of meteorology andsnow pit observations at Tohkamachi Experimental Station, Forestry andForest Products Research Institute, Japan (IX) (2014–15 to 2018-19, five
winter periods). Bulletin of Forestry and Forest Products Research Institute18(4), 393–443. doi: 10.20756/ffpri.18.4_393, (in Japanese).
Takeuchi Y, Nishimura K and Patra A (2018) Observations and numericalsimulations of the braking effect of forests on large-scale avalanches.
Annals of Glaciology 59(77), 50–58. doi: 10.1017/aog.2018.22

14
Techel F, Müller K and Schweizer J (2020) On the importance of snowpackstability, the frequency distribution of snowpack stability, and avalanche sizein assessing the avalanche danger level. The Cryosphere 14, 3503–3521. doi:
10.5194/tc-14-3503-2020
Teich M, Bartelt P, Grět-Regamey A and Bebi P (2012b) Snow avalanches inforested terrain: influence of forest parameters, topography, and avalanchecharacteristics on runout distance. Arctic, Antarctic, and Alpine Research 44(4), 509–519. doi: 10.1657/1938-4246-44.4.509
Teich M, Marty C, Gollut C, Grêt-Regamey A and Bebi P (2012a) Snow
and weather conditions associated with avalanche releases in forests:
rare situations with decreasing trends during the last 41 years. Coldhttps://doi.org/10.1017/jog.2022.85 Published online by Cambridge University Press

Yuta Katsuyama and others
Regions Science and Technology 83–84, 77–88. doi: 10.1016/j.coldregions.
2012.06.007
Weir P (2002) Snow Avalanche Management in Forested Terrain. Victoria, BC:
Research Branch, B.C. Ministry of Forests. Available at https://www.for.gov.
bc.ca/hfd/pubs/Docs/Lmh/Lmh55.htm.
Yamaguchi S, Sato A and Lehning M (2004) Application of the numericalsnowpack model (SNOWPACK) to the wet-snow region in Japan. Annals
of Glaciology 38, 266–272.
Zeidler A (2004) Forecasting skier-triggered avalanches in the Columbia mountains of Canada (Ph.D. thesis). Department of Civil Engineering, Universityof Calgary, Calgary, Alberta, Canada.