Probability assessment of slope instability in seasonally cold regions under climate change

To semi-quantitatively assess the effects of climate change on the slope failures, this paper proposes an effective approach for evaluating the influences of climate change on slope stability in seasonally cold regions. To discuss climate change, this study firstly analyzes the trend of the two main climate factors (precipitation and air temperature) based on the regression analysis results of the meteorological monitoring data during the past 120 years in different scales (e.g., world, country (Japan), and city (Sapporo)), and the downscaled outputs of three different regional atmospheric models (RAMs) with lateral boundary conditions from three different general circulation models (GCMs). Next, to discuss the effects of different climate factors (air temperature, precipitation, etc.) and to determine the key climate factors on the slope instability, a slope stability assessment approach for evaluating the effects of climate changes on slope instability is proposed through the water content simulation and slope stability analysis with considering freeze-thaw action. Finally, to check the effectiveness of the above assessment approach, assessment of instability of an actual highway embankment slope with the local layer geometry is done by applying the past and predicted future climate data. The results indicate that affected by global warming, the air temperature rise in some cold cities is more serious. The climate changes (especially the increase in precipitation) in the future will increase the infiltration during the Spring season. It will lengthen the time that the highway slope is in an unstable state due to high volumetric water content, causing the occurrence of slope failures will be more concentrated in April. While during the Ssummer-Autumn period, the time domain of its occurrence will become wider.


Introduction
In the coldest month of the year, areas where the ground temperature is usually below 0°C, are defined as cold regions [1,2], such as some high latitude regions, like Hokkaido, Japan. In the cold regions, the climate factors (e.g., precipitation and air temperature) and their variations are considered as the significant causes of slope failures or landslides [3,4]. Therefore, many studies have been done to discuss the failure mechanisms and trigger factors of slope failures in the cold regions [5,6].
Ishikawa et al. [7] described that there is a difference in the failure mechanism of slopes between cold regions and warm-temperate regions due to additional factors, i.e., freeze-thaw action. Accordingly, Ishikawa et al. [8] divided slope failures in cold regions into two main classes: slope failures occur during the snow-melting season (March-April-May) and slope failures occur during heavy rainfall season (August-September). Therefore, it is necessary to investigate the influences of climate change on the slope stability in seasonally cold regions.
During the heavy rainfall season, the frequency of severe rainfall events is expected to increase in many mountain areas as a consequence of the influence of global warming [9]. The rainfall infiltration caused decrease in matrix suction is considered as a main cause of slope failures or landslides [10][11][12]. For example, according to the statistical data surveyed by Japan Meteorological Agency [13] for 1000 locations in Japan, the annual average number of rainfall intensity greater than 50 mm/h increased from 195 times (1987~1998) to 226 times (1999~2010). On the other hand, with the intensification of global warming in the past few decades, during the snow-melting season, the occurrences of slope failures in seasonally cold regions are more frequently reported. For example, surveyed by the Ministry of Land, Infrastructure, Transport, and Tourism (MLIT), in April, the average number of slope failures in Japan increased from about 50 times (2008~2012) to near 100 times (2013~2017), almost doubled. While in August, this number increased from about 120 times (2008~2012) to near 240 times (2013~2017), also doubled. Therefore, the influences of climate change on slope stability are intensifying. However, little attention is paid on the influences of climate changes on slope failures in cold regions [11,14,15].
From the view of climate changes, many attempts were also made to assess the influences of climate change on slope failures by exploiting downscaled synthetic future climate scenarios obtained from the GCMs as an input of slope stability [9,16,17]. Chiang and Chang [9] pointed out that the average annual maximum rainfall will have a 15.2% increase in 2010-2099 compared with that in 1960-2008, and the average total unstable area will have a 12.8% increase in Taiwan, China. Ciabatta et al. [17] employed the PRESSCA system (an operational early warning system) for evaluating the influences of climate change on landslide hazard in the Umbria territory, Italy for three different periods (1990-2013 (baseline), 2040-2069, and 2070-2099). They found that in the future period, mainly during the winter season, a general increase (up to more than 40%) in landslides occurrence is expected for the Umbria territory, Italy.
Accordingly, the objectives of this study are (1) to semi-quantitatively assess the effects of climate changes on the uncertainty of climate factors, and (2) to discuss the influences of climate changes on the occurrence of slope failure and its probability distribution. For these purposes, firstly, this study analyzes the trend of the two main climate factors (precipitation and air temperature) in the past 120 years and in the future based on the regression analysis and hybrid-downscaling method. Next, the effects of climate changes on slope instability are semi-quantitatively assessed by using a proposed slope stability assessment approach through the water content simulation and slope stability analysis using a 2dimensional (2D) FEM homogeneous slope model with considering freeze-thaw action. Finally, to check the effectiveness of the above assessment approach, assessment of instability of an actual highway embankment slope with the local layer geometry is done by applying the past and predicted future climate data.

Slope stability assessment approach
Modelling freeze-thaw action and long-term moisture fluctuation Slope failures in cold regions are usually subjected to freeze-thaw actions and rainfall and/or snowmelt water infiltration. Siva Subramanian et al. [18] proposed a slope stability assessment approach to simulate the longterm fluctuation of the moisture and factor of safety (FOS) for cold region slopes based on non-isothermal seepage flow and limit equilibrium method (LEM) as shown in Fig. 1. In this assessment approach, the influences of freeze-thaw action and long-term moisture fluctuation on the slope stability are taken into consideration. Moreover, the rainfall and/or snowmelt water infiltration can be simulated. However, runoff is ignored as this assessment approach is established by using the 2D GeoStudio software [19].

Modelling atmosphere-ground interactions
In order to simulate the influences of climate change on slope stability, the interactions between atmosphere and ground (heat exchange, water transfer, etc.) need to be modeled. The mass balance equation could be used for calculating the water flux on the ground surface [19]: where, P indicates precipitation (mm/h); SM indicates snowmelt water (mm/h); I indicates infiltration (mm/h); R indicates runoff (mm/h) and α is the slope angle (°). Snow accumulation models often use average air temperature to distinguish the precipitation as rainfall or snowfall [20,21]. The threshold value of air temperature to differentiate rainfall and snowfall is 0°C and snowmelt is assumed to occur when snowpack is present on the ground with the positive air temperature. The snowfall is equal to the precipitation under negative air temperature [19].
where, SF is snowfall (mm/h); f s is snow fraction, and SCF is snow correction factor with the value of 1 [22].
The snow fraction, f s is defined as, where, T a is average air temperature (°C) and T f is freezing point temperature (0°C). When air temperature is negative, snow will accumulate on the surface of the soil. Once the air temperature is positive, the snow will melt into water then penetrate into soil. The amount of melted water is calculated in the form of the snow water equivalent.
where, ΔSWE is difference in snow water equivalent (mm/h).

Governing equations for non-isothermal seepage flow
The 2D seepage flow is governed by the equation given by Richards [23] and Kim and Heydinger [24].
where, k x (ψ) and k y (ψ) indicate unsaturated hydraulic conductivities in x-direction and y-direction (m/s), respectively; ψ indicates pore water pressure head (m); c(ψ) indicates specific moisture capacity (1/m). The 2D thermal flow is governed by the equation given by Kim and Heydinger [24].
where, T is temperature (°C); λ T is thermal conductivity (kJ/(Day•m•°C)); v x and v y are the velocities of the water flow in x and y directions (m/s); C w is the heat capacity of water (kJ/(m 3 •°C)); C is the volumetric heat capacity (kJ/(m 3 •°C)).

Governing equation for unsaturated soil shear strength
The shear strength of an unsaturated is expressed based on Bishop's effective stress principle by Vanapalli et al. [25] as given by.
where, τ is shear strength of unsaturated soil (kPa), σ n is net total stress (kPa), u a is pore air pressure (kPa), u w is pore water pressure (kPa), c' is effective cohesion (kPa), ϕ' is effective angle of internal friction (°) and χ is the matrix suction coefficient which varies from 0 to 1 depending on the volumetric water content (θ w ), saturated volumetric water content (θs) and the residual volumetric water content (θr).
Meteorological monitoring data analysis and future climate prediction Meteorological monitoring data analysis To semi-quantitatively assess the influences of climate changes on the uncertainty of climate factors in different scales (e.g. world, country (Japan), and city (Sapporo, a city located in the seasonally cold region of northern Japan)), this study analyzes the trend of the annual average air temperature deviation ( Fig. 2(a)), annual cumulative precipitation deviation ( Fig. 2(b)), magnification of precipitation ( Fig. 2(c)), and the number of days with air temperature below zero ( Fig. 2(d)) by the regression analysis of the meteorological monitoring data of the past 120 years (1900-2020) of the world, Japan, and Sapporo referred from JMA. The benchmark of deviation is the average value of 30 years (1981-2010). The annually average air temperature deviation shown in Fig. 2(a) and annually cumulative precipitation deviation shown in Fig. 2(b) are the difference between them of each year and this benchmark. Figure 2(a) shows that the world annually average air temperature increases approximately 0.8°C per 100 years, while the values of Japan and Sapporo are 1.27°C per 100 years and 2.46°C per 100 years, respectively. It means that global warming has a significant impact on Japan, especially on high latitude area or seasonal snowy cold regions in Japan like Sapporo. Figure 2(b) plots the annually cumulative precipitation deviation. Six intervals of 20 years are divided over the past 120 years to analyze the change of the annual cumulative precipitation deviation. It can be seen that the number of positive values on the three scales are 12, 14, and 15, respectively, in the past 20 years (2000-2020), which are all the maximum values since 1920. It means from 1920, the number of years with heavy rainfall is increasing. Figure 2(c) and Fig. 2(d) shows the magnification of precipitation (the ratio of annually cumulative precipitation of each year to the benchmark) and the number of days with air temperature below zero in Sapporo, respectively. The results of fitting curves display that the magnification of precipitation increases by 0.094 per 100 years, and the number of days with air temperature below zero decreases 16.8 days per 100 years in Sapporo. It means that the precipitation is increasing and winter is getting shorter in Sapporo.
To investigate whether climate change has different effects on air temperature in different seasons, the air temperature change in four seasons, i.e., December-January-February (DJF), March-April-May (MAM), June-July-August (JJA), and September-October-November (SON), in the world, Japan, and Sapporo are analyzed separately as shown in Fig. 3. It is recognized that in all seasons, the influences of climate change on the air temperature of Sapporo are very significant, followed by Japan, compared with the world. It is noted that in Sapporo, climate change has the most obvious effect on the air temperature in the Winter season (DJF, increases by 2.89°C per 100 years), followed by Spring season (MAM, increases by 2.63°C per 100 years), and then Autumn season (SON, increases by 2.30°C per 100 years), and finally Summer season (JJA, increases by 2.03°C per 100 years). The increase in average air temperature in Autumn season and Winter season will cause the thickness of frozen soil to decrease, and in Spring season will cause the increase of snowmelt water, which can be considered as the main reason of the increase in the slope failures in the cold regions in Japan.

Future climate prediction based on dynamical downscaling techniques
To further analyze the future climate change, the future climate data (2090-2091), i.e., meteorological simulation data at Nakayama Pass, Sapporo, Hokkaido, Japan is estimated by using another dynamical downscaling technique, the hybrid-downscaling method proposed by Yamada et al. [26]. The meteorological simulation data is obtained by downscaling the outputs of 3 different RAMs with lateral boundary conditions from 3 different GCMs. Figure 4 shows past meteorological monitoring data (2012-2013) and future meteorological simulation data (2090-2091) at Nakayama Pass. It shows that the average maximum air temperature increases 0.6°C, and the average minimum air temperature increases by 3.29°C in 2090-2091. Therefore, the average air temperature increases by about 1.945°C in the next 80 years at Nakayama Pass, which is consistent with an increase of 2.46°C per 100 years in Sapporo obtained from the regression analysis results of the meteorological monitoring data ( Fig. 2(a)). Besides, the annual cumulative precipitation increases from 1882 mm (2012-2013) to 2131 mm (2090-2091). The magnitude of the precipitation increases by about 0.13, which is slightly larger than the value of the meteorological analysis data (0.094 per 100 years) in Sapporo shown in Fig. 2(c).

Instability assessment of homogeneous conceptual slope under climate change
According to meteorological analysis data obtained from the regression analysis and the meteorological simulation data obtained from hybrid-downscaling method, climate change is gradually increasing the precipitation and air temperature in Japan (especially in cold regions).
To assess the influences of climate change on the occurrence of slope failure and its probability distribution, water content simulation and slope stability analysis were performed using a 2D homogeneous conceptual slope model with the consideration of freeze-thaw actions using the slope stability assessment approach. The applicability and reliability of the slope stability assessment approach in the analysis of slope stability during the snow-melting season and heavy rainfall season in cold regions has been clearly validated by Siva Subramanian et al. [18].
Numerical model and boundary conditions of 2D numerical homogeneous slope model A 2D numerical model is established to simulate slope stability under climate change. Figure 5 shows the numerical model used in this study and the applied boundary conditions. The numerical models are simulated with different slope heights (5 m, 10 m, 15 m, and 20 m) and different slope angles (30°, 35°, 40°, and 45°) to reflect various in situ conditions of homogeneous conceptual slopes. During the simulation, climate data is applied on the soil surface (ab, bc, cd in Fig. 5). Adiabatic boundary conditions are set at the right side (de in Fig. 5), the bottom side (ef in Fig. 5), and the left side (af in Fig. 5). Due to the left side of the model is assumed to adjacent to the bedrock, no flow boundary is applied on the left side (af in Fig. 5) and free drainage boundary is applied on the right side (de in Fig. 5) and bottom side (ef in Fig. 5). Figure 6 shows the original climate data at Nakayama Pass during 2012-2013, which were obtained from the AMeDAS. According to the regression analysis results of the historical meteorological monitoring data of Sapporo, the assumption for precipitation is 1.2 times increased and for air temperature is 2°C increased in the 80 years later. Therefore, 4 kinds of assumed climate boundary conditions as listed in Table 1 are applied to 16 kinds of numerical models (4 slope heights × 4 slope angles). A total of 64 numerical simulations are conducted for Nakayama Pass. Table 2 lists the soil material properties. Dry density (ρ d ), porosity (n), effective cohesion (c'), effective internal friction angle (ϕ'), and saturated hydraulic conductivity (k s ), were measured by Hokkaido Regional Development Bureau [27]. The SWCC and unfrozen volumetric water content function of soil at Nakayama Pass is plotted in Fig. 7. As shown in Fig. 7(b), the unfrozen volumetric water content is related to only temperature. In other words, the condensation and dissipation of ice in the soil over time and moisture migration are not considered, which is one of the limitations of this study. Therefore, the influence of the time effect and moisture migration on unfrozen volumetric water content in the soil still needs to be further investigated. The soil water characteristic curve (SWCC) of the slope is estimated from the grain size distribution curve by referring to Fredlund et al. [28]. The rest parameters are estimated based on the grain size curve of soil [29]. The unfrozen water content of the soil is estimated using the relationships given by Black and Tice [30].
where, θ uw (T) is unfrozen volumetric water content at temperature T, θ is function of volumetric water content against matric suction as, C f is parameter ranges between 1 to 2.2 for fine to coarse grained soils.

Simulation results of 2D numerical homogeneous slope model
To assess the influences of climate change on the occurrence of slope failures, water content simulation and  Fig. 8. Figure 8(a) displays the calculated and measured snow depth. Seeing Fig. 8(a), the results suggest that the calculated snow depth under case 1 is consistent with the measurement data, meaning that the simulation results are reliable. It is recognized that when the air temperature is increased by 2°C, the time when the snow disappears is advanced from June 1st (case 1) to May 23rd (case 2), and the maximum snow depth is decreased from 2.79 m (case 1) to 2.66 m (case 2). When the precipitation is increased by 1.2 times, the maximum snow depth is increased from 2.79 m (case 1) to 3.42 m (case 3). Figure 8(b) and Fig. 8(c) show that the increase in the air temperature increases snowmelt and infiltration during the Spring season and the end of the Autumn season. The main reason is that on one hand air temperature controls snowmelt during those periods and the snowmelt water is the main water supply of the infiltration. On the other hand, when the average air temperature is increased by 2°C, during Spring season and Autumn season, some negative air temperatures become positive, which causes some snowfalls are converted into the rainfalls as shown in Fig. 8(d). At the same time, more snow on the ground surface melt into water and infiltrate into the soil reducing the slope stability.
To assess the influences of climate change on slope stability subjected to snowmelt and heavy rainfall respectively, a year is divided into two different intervals according to the different causes of slope failures: Winter-Spring period (the key factor of slope failure is snowmelt water) and Summer-Autumn period (the key factor of slope failure is heavy rainfall). The probability distribution is described by probability density function, which is a function that can be used to describe the probability of the slope failures in a time interval (1 month in this study).
where, f(x) is the normal density at x, which is parameterized in terms of the mean and the variance, denoted by σ and μ 2 , respectively. Figure 9 (a) and Fig. 9 (b) show the total number of days with slope failure (FOS < 1) and the fitting normal  density curve (NDC) of numerical results during the Winter-Spring period and the Summer-Autumn period, respectively. Figure 9 suggests that when air temperature increases 2°C (case 2) compared with the original climate conditions (case 1), during the Winter-Spring period, the μ will decrease and the σ will increase. The NDC becomes chunky and moves to March. It means that the increased air temperature causes more slope failures to occur earlier and broadens the time domain of slope failures in the Spring season. The main reason is that the increased air temperature leads to the shortening of Winter season and the advancement of the snow-melting season. However, during the Summer-Autumn period, it is the exact opposite that the increased air temperature increases the value of μ and decreases the value of σ. The NDC becomes thin and tall, and moves to September. It means that the increased air temperature causes more slope failures to occur later and narrows the time domain of slope failures in the Autumn season. When precipitation increases 1.2 times (case 3) compared with the original climate conditions (case 1), during the Winter-Spring period, the μ will increase and the σ will decrease, and during the Summer-Autumn period, μ will decrease and the σ will increase. It means that the increased precipitation causes more slope failures to occur later and narrows the time domain of slope failures in the Spring season, and more slope failures to occur earlier and broadens the time domain of slope failures in the Summer season.
When considering changes in air temperature and precipitation together, it can be seen that during the Winter-Spring period both μ and σ decrease, while during the Summer-Autumn period both μ and σ increase. This is consistent with the trend of slope failures in Japan investigated by the MLIT from 2008 to 2017 as shown in Fig. 10. It indicates that during the Winter-Spring period the occurrence of slope failures will be more concentrated in April, while during the Summer-Autumn period, the time domain of its occurrence will become wider. Therefore, the numerical results are believable through the verification of the investigation results. It indicates that the meteorological analysis data are reliable and the proposed assessment approach is available for evaluating the effects of climate changes on slope instability.

Assessment of actual highway slope instability under climate change
In previous chapter, an assessment approach for evaluating the effects of climate changes on slope instability is proposed using the 2D homogeneous conceptual slope model with considering different slope heights and different slope angles. To check the effectiveness of the above assessment approach, assessment of instability of an actual highway embankment slope with the local layer geometry is done by applying the past (meteorological monitoring data during 2012-2013) and predicted future climate data by dynamical downscaling techniques (meteorological simulation data during 2090-2091). Along the national highway route 230 in Hokkaido, Japan, at Nakayama Pass, the slope failure occurred successively on May 4th, 2012, and April 7th, 2013 as shown in Fig. 11. On May 4th, 2012, within the range of about 40 m in length and up to 110 m in transverse width, approximately 13,000 m 3 of soil collapsed and caused 20 days road closure. On April 7th, 2013, another slope failure occurred with more than 11,000 m 3 of sediment collapsed [31]. The increase in the snowmelt water due to climate change was considered as the main causes of the slope failures.

Numerical model and boundary conditions of the highway slope at Nakayama pass
The 2D numerical model with the slope stratigraphy for highway slope at Nakayama Pass is established, which consists of three layers of embankment filling, talus, and bedrock as shown in Fig. 12. The parametric analysis of slope failures at Nakayama Pass is performed by inputting meteorological monitoring data during 2012-2013 (measured by AMeDAS) and meteorological simulation data during 2090-2091 (simulated by GCMs) as shown in Fig. 4. The measured groundwater table is given in Fig. 12. Table 3 lists the five kinds of boundary   conditions for parametric studies. The soil properties are given in Table 4. The bedrock is modeled as an impenetrable material.
Simulation results of the highway slope at Nakayama pass To check the effectiveness of the above assessment approach, assessment of instability of an actual highway embankment slope with the local layer geometry is done by applying the past and predicted future climate data. The simulation results of the variation in the FOS, ground surface temperature, volumetric (vol.) water content and infiltration are given in Fig. 13. Figure 13(a) displays that during snow-melting season, the FOS has a rapidly decrease due to the snowmelt water infiltration increases the vol. water content as shown in Fig. 13(d), which causes the slope to be in an unstable state (FOS < 1). After the snow-melting season, the FOS gradually rises to the normal level. It is noted that there is a much longer (more than doubled) unstable state of the highway slope at Nakayama Pass under B3 (53 days) and B5 (64 days) boundary conditions compared with the B1 (22 days) boundary condition. However, the unstable state of the highway slope at Nakayama Pass has almost no significant increase under B2 (23 days) and B4 (25 days) boundary conditions. It means that climate change (especially the increase in precipitation) in the future will enhance the instability of the slopes and lead to more slope failures occur. Figure 13 under B3 and B5 boundary conditions with an increase of approximately 69% and 80%, respectively. During this period, the ground surface is unfrozen as shown in Fig. 13(c). Therefore, the infiltration increases the vol. water content eventually causes slope failure (see Fig. 13(d)). Accordingly, it can be obtained that climate change (especially the increase in precipitation) in the future will increase the infiltration and cause more slope failures during the Spring season.

Discussions
To assess the effects of climate change on the slope stability, this paper proposes an assessment approach for evaluating the effects of different climate factors (precipitation, air temperature, wind speed, etc.) on slope instability through the water content simulation and slope stability analysis with considering freeze-thaw action in seasonally cold regions. The proposed assessment approach can evaluate the effects of the climate change on the slope stability by analyzing the change of the shape and translation of the normal density curve (NDC) of the distribution of slope failures in 1 year, which provides an effective method to semi-quantitatively evaluate the impact of climate change on the occurrence time and distribution of slope failures. During the analysis, three kinds of different climate data are discussed: meteorological monitoring data, meteorological analysis data, and meteorological simulation data. The meteorological monitoring data is the direct monitoring data so it is more reliable. However, this kind of climate data is historical data so it can not be used in the future. The meteorological analysis data is the assumed future climate data that add the average changes in precipitation and air temperature from the regression analysis to the meteorological monitoring data. Therefore, this kind of climate data can be obtained very easily and quickly, while the differences of climate change in different seasons are not considered. The meteorological simulation data are obtained from GCMs. The differences of climate change in different seasons are considered, while this kind of data is more complex and difficult to obtain. In addition, the analysis of the slope stability in this study is based on a 2D FEM model, therefore, extending the simulation model to 3D is the future assignment of this study.

Conclusions
According to the regression analysis results of the meteorological monitoring data during the past 120 years in different scales (e.g., world, country (Japan), and city (Sapporo)), affected by global warming, the air temperature rise in some cold cities is more serious. Meanwhile, the precipitation is also increasing and Winter is getting shorter in these cold regions. In different seasons, climate change has the more obvious influences on the increase in the air temperature in the Winter season (DJF) and Spring season (MAM)   Climate change (especially the increase in precipitation) in the future will increase the infiltration during the Spring season. It will lengthen the time that the highway slope is in an unstable state, thereby enhancing the instability of the slopes and lead to more slope failures occur in the future. The increased air temperature leads to the shortening of winter season and the advancement of the snow-melting season. It causes the slope failures occur earlier and broadens the time domain of slope failures in the Spring season, and causes more slope failures to occur later and narrows the time domain of slope failures in the Autumn season. The increased precipitation causes more slope failures to occur later and narrows the time domain of slope failures in the Spring season, and causes more slope failures to occur earlier and broadens the time domain of slope failures in the Summer season. When considering changes in precipitation and air temperature together, during the Winter-Spring period, the occurrence of slope failures will be more concentrated in April, while during the Summer-Autumn period, the time domain of its occurrence will become wider.
The research findings propose an effective method to semi-quantitatively evaluate the influences of climate change on slope stability and the probability distribution of the slope failures by using numerical simulation, regression analysis, and statistical analysis. It can be used in the adjustment of counter measures of slope failures at the actual design and maintenance works in the future.