To study the influencing factors of the distribution and variation trend of snow BC concentration in northeast China, the spatial distribution data of China's elevation (DEM) released by the resource and environment data cloud official website (http://www.resdc.cn/ data.aspx?DATAID=284) (accessed on 10 May 2019), the average wind-speed data released by the scientific data center in cold and arid regions (http://westdc.westgis.ac.cn/data/ 7a35329c-c53f-4267-aa07-e0037d913a21) (accessed on 10 May 2019), and the night light time data of the DMSP satellite (https://ngdc.noaa.gov/eog/dmsp/downloadV4composites. html) (accessed on 12 May 2019) were used to discuss the effects of topographic factors, population density, energy factors, and meteorological conditions on the distribution and variation trend of snow BC concentration in northeast China.

#### 2.2.2. Field Data

To validate the model, twice field samples were taken from January to February 2010 and January 2014, totaling 28 points, including 19 in 2010 and 9 in 2014 (Figure 2). The specific sampling process was the following: clean stainless-steel utensils were used to dig out snow pits and make a vertical profile in one direction. When sampling, according to the order of "left" and "right", two snow samples of 0–5 cm cover in vertical section were collected and placed in Whirlpak sampling bags, sealed, and put into an incubator to ensure that they remained in a solid-state until the next processing step in the temporary laboratory. In the temporary laboratory, samples of the collected snow were scooped into clean glass beakers with clean stainless-steel containers and quickly melted in a microwave oven. After measuring the volume of the melted samples, the melted snow water was filtered through a 0.4-mm Nuclepore filter by a manual pump immediately to separate BC and other water-insoluble substances [46], and then an ISSW spectrophotometer was used to analyze the BC and other light-absorbing substances in the filter samples [25,47]. Since the BC concentration was estimated based on the Ångström absorption index of the sample on the filter, the BC concentration of the sample given was a range, not a definite value [48]. For example, the uncertainty range of the BC concentration of a snow sample in Changchun City in early February 2010 was 250–650 ng/g.

**Figure 2.** Locations of the 28 sampling sites (blue is the sampling point in 2010 and green is the sampling point in 2014), ranging in latitude from 42◦ to 51◦ N and in longitude from 120◦ to 131◦ E.

#### *2.3. Methods*

#### 2.3.1. Snow Recognition Method

Snow pixel judgment was the basic data for simulating snow BC concentration. In this study, we used the SNOWMAP algorithm to identify snow-covered pixels, which was proposed by Hall et al. [49], and subsequently improved for the production of MODIS snow cover products [50] and has been validated and applied worldwide [51,52].

The core concentration of the SNOWMAP algorithm is normalized difference snow index (NDSI), which uses the spectrum characteristic of the high reflectivity of snow in the visible light band and low reflectance in the near-infrared band to distinguish between clouds and snow. The calculation method is

$$NDSI = \frac{R(b4) - R(b6)}{R(b4) + R(b6)} \tag{1}$$

where, *R*(*b*4) and *R*(*b*6) representrepresents the surface reflectance of band 4 (545–565 nm) and band 6 (1628–1652 nm) in the MOD09GA data, respectively.

Hao et al. [53] found that the NDSI threshold was 0.4. The NDSI of water is also greater than 0.4, so we eliminated the influence of water on snow pixel based on the characteristic that the reflectivity of snow in the second waveband was greater than, or equal to, 0.11, and the reflectivity of water in this waveband was not greater than 0.11 [50]. Therefore, the discriminant formulas of snow pixels in this study were

$$\begin{cases} \begin{array}{c} NDSI \geq 0.4\\ R(b2) > 0.1 \end{array} \end{cases} \tag{2}$$

where, *R*(*b*2) represents the surface reflectance of band 2 (841–876 nm).

In this paper, the snow pixels identified by the SNOWMAP algorithm were marked as 1, and the other non-snow pixels were marked as 0. A snow cover map with a spatial resolution of 500 m × 500 m was generated.

#### 2.3.2. ART Model

The ART model is a simplified radiative transfer model developed by Kokhanovsky et al. [54,55]. The basic principle is to regard snow particles as irregular particles, consider the effect of pollutants on the absorption characteristics of the snow surface, the optical characteristics of a single particle are calculated using the principles of geometrical optics, and then the asymptotic analysis solution of the radiative transfer model is obtained using the asymptotic analysis method. The model driving parameters mainly include visiblelight reflectance, near-infrared reflectance, solar zenith angle, viewing zenith angle, and relative azimuth angle. Snow characteristics, such as the grain effective size (*αef*) and BC concentration (*c*), can be estimated on regional and global scales.

The formula for simulating the BC concentration of each grid in snow based on the ART model is:

$$R\_1 = R\_0 \exp\left[-\frac{4f(\mu\_\prime \mu\_0, q)}{\sqrt{3(1 - g)}} \sqrt{\frac{2}{3} B a\_{s, 1} c a\_{cf}}\right] \tag{3}$$

$$R\_2 = R\_0 \exp\left[-\frac{4f(\mu\_\prime \mu\_{0\prime} \,\, q)}{\sqrt{3(1-g)}} \sqrt{a\_{i,2} \beta \mathcal{K} \alpha\_{ef} + \frac{2}{3} B \alpha\_{s,2} c \alpha\_{ef}}\right] \tag{4}$$

There are 15 parameters in the equations (Equation (3) and (4)), and the physical meaning and units are listed in Table 1. The snow BC concentration in this study is relative snow concentration.

**Table 1.** Parameters used in ART model.

*R*1, *R*2, *ϑ*0, and *ϑ* in equations (Equations (3) and (4)) can be downloaded from the MODIS data. Other parameters can be obtained using the following equations.

1. Reflection function of the semi-infinite snow layer *R***0**

*R*0 was calculated by the following equation (Equation (2)) [56,57]

$$R\_0(\mu, \mu\_0, \varphi) = \frac{A + B(\mu + \mu\_0) + \mathbb{C}\mu\mu\_0 + P(\theta)}{4(\mu + \mu\_0)}\tag{5}$$

The parameters in Equation (5) are listed in Table 2.

**Table 2.** Parameters needed to calculate *R*0.

The phase function *P*(*θ*) does not depend on the wavelength, therefore the function *R*0 (*μ*, *μ*0, *φ*) is the same in the entire spectrum. The phase function *P*(*θ*) is calculated as:

$$P(\theta) = 11.1 \exp(-0.087\theta) + 1.1 \exp(-0.014\theta) \tag{6}$$

where, *θ* is defined as *θ* = *a* cos(−*μμ*0 + *ss*0 cos *ϕ*), *s* = sin(*ϑ*), *s*0 = sin(*ϑ*0), *μ*0 = cos(*ϑ*0), and *μ* = cos(*ϑ*).

#### 2. Parameter *f*(*μ*, *μ*0, *φ*)

*f*(*μ*, *μ*0, *φ*) is determined by the escape function and the semi-infinite space snow reflectance function, and its calculation formula is as follows:

$$f(\mu, \mu\_0, \boldsymbol{\varrho}) = \frac{\boldsymbol{\mu}(\mu\_0)\boldsymbol{\mu}(\boldsymbol{\mu})}{R\_0(\mu, \mu\_0, \boldsymbol{\varrho})} \tag{7}$$

where, *u*(*μ*0) is called the escape function, which is solved according to the following empirical formulas:

$$
\mu(\mu\_0) = \frac{3}{7}(1 + 2\mu\_0), \; \mu(\mu) = \frac{3}{7}(1 + 2\mu) \tag{8}
$$

#### 3. Absorption coefficient of ice crystal and BC

In the visible light band (channel 1), the absorption of light by the snow layer is mainly caused by pollutants (mainly BC); in the near-infrared band (channel 2), the absorption of light by the snow layer is mainly caused by ice crystals and BC.

The formula of the absorption coefficient of BC is

$$\mathfrak{a}\_{s,i} = 4 \prod \chi\_s(\lambda) / \lambda\_i \tag{9}$$

where, *i* = 1,2, indicating the first and second channels. Since the refractive index of BC is independent of wavelength, *χs*(*λ*) is a constant of 0.46. *λ* is the wavelength; that is, *λ*1 = 460 nm and *λ*2 = 865 nm.

The expression of the light absorption coefficient of ice crystals is

$$
\pi\_{\rm icc,2} = 4 \prod \chi\_{\rm icc}(\lambda) / \lambda \tag{10}
$$

where, *αice*,2 is the ice absorption coefficient at the wavelength *λ*2; *χice*(*λ*) is the imaginary part of the refractive index of ice, which is a function of wavelength. Here, we use *χice*(*λ*) = 2.5 × 10−9.

4. Relative azimuth angle *φ*

*φ* is calculated with Equation (11):

$$
\varphi = saa - vaa
\tag{11}
$$

where *saa* is the solar azimuth angle and *vaa* is the viewing azimuth angle.

Entering all the above parameters into the ART model, we simulated the BC concentration of each snow pixel with a spatial resolution of 500 × 500 m.

Based on the above data and formula, we analyzed the BC concentration of snow cover in northeast China, and the research framework is shown as follows (Figure 3):

**Figure 3.** The framework of the study.

#### **3. Results**

#### *3.1. Model Verification*

We verified the accuracy of the simulated values by comparing the simulated snow BC concentration with the measured values. Since the measured value was an uncertain range, as long as the simulated value was distributed within the range of the measured value we could consider the simulation result credible. We also used the literature comparison method to verify the simulated value. The field sampling measured datum were used to verify the simulation values of the same period (Figure 4). There are significant disagreements in the measured and modeled snow BC concentrations at individual sites, but over the region as a whole there is no apparent systematic bias in the modeled vs. measured concentrations. Zhao et al. [58] simulated the BC concentration of snow cover in northeast China in 2010, and the simulated value was above 1000 ng/g, similar to the results of this study (1197.468 ng/g). For areas north of 48◦ N, the simulated value was higher than the measured value. However, Flanner et al. [36] simulated the snow BC content in the high latitudes of northeast China, and the simulation result was above 1000 ng/g, which is similar to the simulation result of this study. Based on the above comparison, the ART model has a high ability to simulate the BC concentration of snow cover in northeast China, and the simulation results are credible.

#### *3.2. Annual BC Concentration*

#### 3.2.1. Temporal Variation

From 2001 to 2016, the average annual BC concentration of snow cover in northeast China was between 1098.927 and 1257.300 ng/g, with an average value of 1197.468 ng/g. The maximum value appeared in 2016 and the minimum value appeared in 2002 (Figure 5). No significant (*p* > 0.05) increase in BC concentration was observed from 2001 to 2016, with an annual growth rate of 5.137 ng/g.

**Figure 4.** BC concentrations in the top snow layer from the campaign observations (black) and the corresponding ART simulations (red). The column plot of observations shows the minimum and maximum possible values from estimates; the triangle within the column shows the average of observations; and the circle within the column shows the model results.

**Figure 5.** Average annual variability of BC concentration on snow cover in northeast China from 2001 to 2016.

#### 3.2.2. Spatial Distribution and Variation

Figure 6 shows the spatial distribution of the average BC concentration over northeastern China during 2001 to 2016. There was generally a strong spatial difference in BC concentration that varied considerably over the 310–1561 ng/g range in the study area. Generally, there were two high-value (>1200 ng/g) regions of BC concentration in northeast China, both of which were located in Heilongjiang province. One is an industrial corridor composing Harbin, Daqing, Qiqihar, Suihua, and Heihe; and the other is a coal-forest industrial area composing Hegang, Jiamusi, and Shuangyashan. Conversely, the average BC concentration values were largely <850 ng/g in the mountainous areas of the study area, where human activity is minimal.

**Figure 6.** Spatial distributions of the averaged BC concentration over the study area during December to February 2001–2016.

From 2001 to 2016, the BC concentration in 81.805% of the study area showed an increasing trend (Figure 7). Approximately 6.975% of the study area showed a significant (*p* < 0.05) increase, which was mainly distributed in the northeast of the Daxinganling region, east of Heihe, north of Qiqihar, northeast of Harbin, north of Suihua, and southeast of Shuangyashan. In contrast, only 18.195% of the study area showed an insignificant (*p* > 0.05) decrease in BC concentration, which was mainly distributed in the eastern part of Heilongjiang Province and the central and eastern parts of Liaoning Province. Notably, by combining Figures 5 and 6, we found that the heavy snow BC polluted areas centered on the southern part of the Heihe had significantly worsened pollution. In addition, areas where snow cover was less polluted by BC, such as the Xiaoxing'an Mountains, Changbai Mountains, and southern Liaoning, were also becoming more polluted. These results indicate that the BC pollution degree of snow cover in northeast China increased overall from 2001 to 2016, especially in heavily polluted areas.

#### *3.3. Monthly Snow BC Concentration*

#### 3.3.1. Temporal Variation

The average BC concentration was the highest in December (1344.588 ng/g), and it was gradually decreased in January (1248.619 ng) and February (983.635 ng/g), which was related to the cleaning effect of snowfall on the air. The BC concentration in December and January showed an increasing trend from 2001 to 2016, with annual increases of 7.087 ng/g and 6.556 ng/g, respectively (Figure 8); while in February, it showed a decreasing trend, with an annual reduction of 0.079 ng/g. The coefficients of variation of BC concentration from December to February were 5.601, 8.392, and 8.594, respectively, indicating that the inter-annual distribution of BC concentration in each month was relatively stable.

**Figure 7.** Trends of BC concentration in the top snow layer in northeast China from 2001 to 2016.

**Figure 8.** Average annual variability of BC concentration on snow cover in December, January, and February from 2001 to 2016.

#### 3.3.2. Spatial Distribution and Variation

Figure 9 shows the spatial distribution of the average BC concentration in December, January, and February over northeastern China during 2001 to 2016. There was generally a strong spatial difference in BC concentration monthly that showed higher values in the north than in the south. In December, the BC content was greater than 850 ng/g, and there were three high-value centers (>1400 ng/g) in the spatial distribution, located in the Daxinganling region, northern Qiqihar and southern Heihe, and Hegang and Jiamusi regions. In January, the high-value center range shrank sharply, especially in northern Qiqihar and southern Heihe. By February, the BC content in the study area was drastically reduced, falling below 1200 ng/g. Areas with low BC concentration (<900 ng/g) were located in Yichun City, Huludao City, Chaoyang City, Fuxin City, and Dalian City.

**Figure 9.** Monthly mean value of BC concentration on snow cover in the northeast of China from 2001 to 2016: (**a**) December, (**b**) January, and (**c**) February.

From 2001 to 2016, the increasing trend of BC content in December and January accounted for 82.519% and 82.789% of the total area, respectively (Figure 10). Among them, the BC concentration in central Heilongjiang Province and southern Liaoning increased significantly (*p* < 0.05), accounting for 9.922% and 10.603% of the total area, respectively, indicating that BC concentration in central Heilongjiang Province increased significantly in December and January. On the contrary, the area with a decreasing trend of BC content in February accounted for 60.110% of study area, and BC concentration decreased significantly (*p* < 0.05) in the southern Greater Xing'an Mountains and northern Heihe, accounting for 2.015% of the entire northeast area.

**Figure 10.** Monthly average snow BC concentration tendency rate of northeast China from 2001 to 2016: (**a**) December, (**b**) January, and (**c**) February.

#### *3.4. Influencing Factors*

#### 3.4.1. Terrain

The BC on the snow is mainly the dry deposition of BC aerosols in the atmosphere, therefore the difference in the spatial distribution of BC concentration on the snow is

mainly determined by the difference in the spatial distribution of the BC aerosols in the atmosphere. Lv et al. [59] found that elevation is the most important factor contributing to the concentration distribution of BC aerosols, followed by population density, and heavy industry scale. Owing to the sparse population of high-altitude areas, the lack of large-scale heavy industrial production, and the high coverage of surface vegetation in these areas, the concentration of BC aerosols is low. The northeast of China is surrounded by mountains on three sides (Figure 11), while the central and eastern regions are plains and terraces. Therefore, the concentration of BC aerosol in the central region is relatively high. This is consistent with the conclusion that the BC content of snow in the central plain area is higher than that in the surrounding mountainous areas. It shows that geomorphology types are one of the factors that affect the BC concentration in snow.

**Figure 11.** Geomorphology types of northeast China.

#### 3.4.2. Energy Consumption

BC aerosols in the atmosphere mainly come from the incomplete combustion of fossil fuels and biomass fuels. Studies showed that population density and the scale of heavy industry are also important factors affecting the distribution of BC aerosol concentration [59]. The value of the night light time index can reflect the amount of power consumption, representing information such as city size, population density, and energy use [60,61]. Therefore, this study selects the spatial distribution of nighttime light index data in northeast China in 2001 and 2013 (Figure 12), which represents information such as population density and energy usage. The high values of night light hours in 2001 and 2013 were located in Shenyang, Changchun, Harbin, and Daqing, showing a clear zonal distribution. The Songnen Plain in the west and the Sanjiang Plain in the east are also areas with relatively high nighttime light index, suggesting that they are the areas with large cities, concentrated industrial production, and high population density in the three northeastern provinces. Comparing the spatial distribution of snow BC in northeast China in 2001 and 2013, it can be seen that the high value area of the nighttime light index in northeast China in 2001 and 2013 is basically consistent with the high value area of snow BC concentration. Comparing the spatial changes of the nighttime light index in 2013 and 2001, the areas with light index in western, central and eastern Heilongjiang Province have increased significantly. This is the main reason why BC concentration of winter snow in Heilongjiang province became the largest region in northeast China from 2001 to 2016, and BC concentration of winter snow in central and eastern Heilongjiang Province increased significantly. In addition, it can also be seen that the density and scale of cities in the southern part of northeast China

have increased significantly. Therefore, the concentration of BC in the snow in the southern part of northeast China in 2013 showed an increasing trend in most areas compared with 2001. The above results indicate that population density, city size, and energy consumption are important factors affecting the concentration of BC in snow.

**Figure 12.** Spatial distributions of Nighttime Lights Time in northeastern China in 2001 and 2013.

#### 3.4.3. Atmospheric Stability

The spatial distribution characteristics of BC concentration on snow cover in northeast China are not only affected by terrain and energy consumption, but also related to winter atmospheric stability. When the atmosphere stratification is in a stable state, the poor diffusion ability of the atmosphere makes it difficult for the pollutants to diffuse; otherwise, the pollutants are easy to diffuse [62]. In this study, the average wind speed was used to characterize the stratification stability of the near-surface atmosphere. The lower the wind speed, the more stable the atmospheric layer structure, and the more likely it is to form high-concentration pollution in local areas. Figure 13 shows the average winter wind speed from 2001 to 2016 (December, January, and February the following year). It can be seen from Figure 12 that the winter wind speed was generally low and the atmospheric stratification was relatively stable. The wind speed in the eastern region was relatively high, mainly due to the formation of canyon winds between the Xiaoxing'anling Mountains and Changbai Mountains. In particular, in the central and northern regions of Heilongjiang Province, the wind speed was the lowest, therefore the atmosphere stratification in this area was stable. In addition, the winter temperature was low and the mixed layer was not high, resulting in frequent occurrence of temperature inversions. These climatic conditions hinder the diffusion of BC aerosols. Moreover, these areas were surrounded by mountains on three sides, and the BC aerosols were more difficult to transport outwards and accumulate near the ground. They enter the snow surface through dry deposition, resulting in northern Qiqihar, southern Heihe, and northern Suihua being the high centers of BC in the northeast.

#### 3.4.4. Land Use Type

Land use type can not only indicate the coverage of surface vegetation, but also a direct result of the impact of human activities on the surface. Comparing the land use types in northeast China in 2020 with the simulation of the spatial distribution of BC in northeast China from 2001 to 2016 (Figure 14), it can be seen that the areas where cultivated land and cities are consistent with high-value areas of BC in the snow, the correlation coefficient between urban area expansion and BC concentration increase is 0.584 (*p* < 0.05). It shows that human activities can also affect the concentration of BC in the snow by changing the land use type.

**Figure 13.** Average winter (December, January, and February the following year) wind speed in northeast China from 2001 to 2016.

**Figure 14.** Spatial distributions of land use type in northeastern China in 2020.

In conclusion, the BC content of snow cover in northeast China is affected by both natural and artificial factors. Natural factors are mainly affected by terrain and stable atmosphere stability. Due to sparse population and no large-scale heavy industry production, the concentration level of BC aerosol in high-altitude areas is relatively low. When the atmosphere stratification is in a stable state, the poor diffusion ability of the atmosphere makes it difficult for the pollutants to diffuse, and the BC concentration increased. Human factors mainly include urban expansion, energy structure change, and land use change, among which urban density, urban size, and population density are important factors affecting the distribution of snow BC concentration. The BC concentration in snow cover is the highest in the region where farmland and city are located, suggesting that human activities affect the spatial distribution of BC content in snow cover by changing the surface type.

#### **4. Discussion**

This study used the ART model combined with remote-sensing data to simulate the BC concentration on snow cover in northeast China. Compared with the existing research results, this study not only analyzed the spatial and temporal distribution of BC concentration on snow cover in northeast China from 2001 to 2016, but also further analyzed the characteristics of snow BC concentration in northeast China during the year and inter-annually, and briefly analyzes its influencing factors.

In this paper, through simulations, a long-term series of snow BC concentrations in a certain area, and the spatial distribution and changes of snow BC concentrations with high spatial resolution year by year were obtained. Compared with previous studies, there are few similar studies, and it is also one of the main contributions of this paper. The research results are of great significance for snow albedo simulation research, improvement of regional climate model simulation accuracy, and snow pollution prevention and control. The results showed that the BC content of snow cover in northeast China did not increase significantly from 2001 to 2016, but increased by 52.299 ng/g from 2012 to 2016, compared with that from 2001 to 2005. The 93.25% of the northeast showed an insignificant trend, and the significant increase areas were mainly distributed in the northeast of the Daxinganling region, the east of Heihe City, the northern edge of Yichun City, and the southeast of Shuangyashan City.

Owing to the influence of local environment, human activities, source emissions, and atmospheric circulation, the global distribution of BC aerosol concentration is highly regional, therefore the BC concentration on snow is also significantly different in space. The results show that in the inaccessible Antarctic, the average BC concentration on snow is only 0.2 ng/g [26], while the BC concentration in the Arctic is higher than that of the Antarctic, with an average of 25 ng/g [24]. At the same time, the results of studies by Clarke and Noone [24], Warren and Clarke [26], Doherty et al. [25], Grenfell et al. [63,64], and McConnell et al. [65] show that the variation of snow BC concentration in the European Arctic are greater than in the Canadian Arctic and Arctic Oceans, and the reason may be related to the distance to the source of BC emissions. The snow BC concentration in the middle and low latitudes of mountainous and urban areas is much higher than that in the Antarctic and Arctic regions. For example, the snow BC concentration in the Alps in December 1992 was 22–302 ng/g [31], and the BC concentration of Olympus in the United States in 2012 was 250 ng/g [66]. The average BC concentration of glaciers in the Tibetan Plateau and Xinjiang from 2001 to 2004 was 63 ng/g [67]. The ice core in the monsoon climate zone of the Tibetan Plateau shows that BC concentration has increased rapidly since the 1980s, and the glaciers in the southeast of the Tibetan Plateau have accelerated to retreat [68], indicating that anthropogenic BC aerosol emissions in southeast Asia play an increasingly important role in glacial melting on the Tibetan Plateau. Mid-latitude urban areas have the highest snow BC concentration. For example, snow samples from southern New Mexico and northern Texas in the United States from 1982 to 1985 show that the BC concentration in urban areas is three times that of remote rural areas [28,29]; the range of snow BC concentration in the city of Detroit from 1984 to 1985 was 17–5700 ng/g [30]. In 2007–2009, the BC concentration in the surface snow reached 2740 ng/g in the Sapporo area of Japan [27]. In January and February 2010 and January 2014, the snow BC concentration of northeast China was 40–1600 ng/g and 50–3700 ng/g, respectively. At the same time, it can be seen from the two sampling results that the snow BC concentration in the northeast heavy industrial area is higher than that in the northern border area [38,48,69]. The above statements all reflect the impact of human activities on regional snow BC concentration.

Since the strong absorption of light by BC aerosols in the atmosphere was discovered in the 1980s, scientists have subsequently discovered that the albedo of visible light in snow depends largely on the concentration of BC on snow. Warren and Wiscombe [3] believe that trace amounts of light-absorbing substances can significantly reduce the snow surface albedo in the visible light band, but have little effect on the albedo in the nearinfrared band. This study summarized previous measurements of simulated snow BC

concentrations and corresponding reduced albedo to obtain a simple equation between the BC concentration and its reduced albedo (Figure 15). Using the equation given in Figure 15, the amount of albedo reduction caused by the BC concentration simulated in this study was roughly calculated (Table 3). The estimation results showed that the larger the snow BC concentration, the more the albedo decreases. These results are only theoretical estimates and require field measurements for future verification.

**Figure 15.** Logarithmic scatter plot for BC concentrations and reduced albedos (hollow circles) in snow and ice based on previous studies, and the simple linear relationship possibly existing between them. The data of the 6 points in the figure are from the literature: 1-Jacobson [8]; 2-Warren [3]; 3-Jacobson [8]; 4-Clarke and Noone [24]; 5-Jacobson [8]; 6-Flanner [7].

**Table 3.** BC concentration and its reduced albedo in the three northeastern provinces and capital cities.

In addition to BC, light-absorbing pollutants on snow cover include organic carbon (OC) and mineral dust (MD), accounting for 11.7% of total light absorption [69]. OC comes from the incomplete combustion of fossil fuels and biomass fuels, which are deposited on the snow surface after dry and wet deposition. MD mainly comes from the arid and semi-arid desert areas of Mongolia. Under the influence of winter monsoon, the dust in the desert area is lifted up and transported to the northeast region and deposited on the surface of the snow there. ART model does not have the ability to detect OC and MD concentration on snow. Therefore, in future research, OC and MD must be taken into account. Based on the optical characteristics of the three substances, BC, OC, and MD, the ART model parameters should be improved to make it more suitable for the characteristics of snow pollution in northeast China.

Compared with the previous research, this study proposed to increase the influence factors of land use type and reached the conclusion that cultivated land and snow BC high value area are highly consistent. But why is the content of BC in the snow in the area where the cultivated land is located so high? The mechanism needs to be further discussed. Because northeast China is one of the three largest black lands in the world, is there a high BC content in the simulation due to incomplete snow cover in the grid and part of the black lands? This also needs to be further revised through field investigations. In addition, this paper used the lighting index to indirectly analyze the impact of population density, energy consumption, and urban expansion on the concentration of BC in the snow. The relationship between the lighting index and population density, energy consumption, and city size and expansion need further analysis. The research conclusions of this study showed that the BC content of northeast China's snow cover increased by 52.299 ng/g from 2012 to 2016 compared with 2001 to 2005. In fact, there was a net loss of population in northeast China in the recent five years. Why is the BC content in snow cover increasing? This also needs to be combined with the actual data for further analysis.

The model only considers the situation that BC is deposited on the snow surface through dry sedimentation, but it does not consider the situation that BC particles dissolve in the snow particles. After BC is deposited in the snow, the snow gradually ages and deforms until the snow melts as time goes by. During this period, the BC particles migrate with the snow melt water, and their effective particle size and optical properties change. In the future, we can observe the change process of BC in snow and perform experimental simulation to establish the model mechanism, so as to more accurately simulate its influence on albedo.

#### **5. Conclusions**

(1) From 2001 to 2016, the average annual BC concentration in northeast China was between 1098.927 and 1257.300 ng/g, with an average value of 1197.468 ng/g. Industrial corridors composed of Harbin, Daqing, Qiqihar, Suihua, and Heihe and coal-forest industrial production areas composed of cities such as Hegang, Jiamusi, and Shuangyashan were high-value centers of snow BC concentration. The area with increased BC concentration accounted for 81.805% of the study area, of which 6.975% was a significant increase area, mainly distributed in the middle of Heilongjiang Province.

(2) The BC concentration of northeast China was the highest in December (1344.488 ng/g), followed by January and February from 2001 to 2016. From 2001 to 2016, the BC concentration of snow cover in December and January in northeast China showed an increasing trend, while the BC concentration in February showed a decreasing trend, with no significant increase or decrease trend in each month.

(3) The high vales of BC concentration in Shenyang, Changchun, Harbin, and Daqing are due to large population density and concentrated industrial production. In the central part of Heilongjiang province, due to the stability of the atmospheric stratum and the influence of topographical conditions, the pollutants were not easy to diffuse, causing Heilongjiang province to become the center of high BC concentration. The increase of cultivated land area, urban density and scale, and the expansion of population area are the main reasons for the increase of BC concentration in snow cover, and also the reasons for the high value of BC content in snow cover in Heilongjiang Province.

This study simulated the BC concentration of surface snow, analyzed the temporal and spatial distribution and variation characteristics of BC concentration, and analyzed the influencing factors of BC concentration of surface snow in northeast China from 2001 to 2016. Compared with previous studies, the main contribution of this study is to analyze the distribution characteristics of in snow cover with large scale and high resolution, and to obtain the variation characteristics of BC concentration in snow (2001–2016) in northeast China. The results of this study provide important basic data for simulating the albedo of snow cover in winter in northeast China and analyzing the influence of snow cover on the regional climate. The shortcomings of current studies lie in the need to further improve the spatial resolution of the simulation of snow BC concentration, refine the analysis of temporal and spatial distribution and variation factors of snow BC, and deepen the study of the influence of BC on snow albedo.

**Author Contributions:** Y.Z. designed the research framework and drafted the manuscript; L.Z. completed the manuscript and made major revision; W.L. was involved in the design and revision of the manuscript; F.Z. and X.Z. were involved in data acquisition and analysis. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China under grant number 41771067 and the Key projects of Natural Science Foundation of Heilongjiang Province under grant number ZD2020D002.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available by request from the corresponding author.

**Acknowledgments:** We sincerely thank Wang Xin of Lanzhou University for providing snowfield sampling data for this study.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

#### *Article* **Spatial Distribution of, and Variations in, Cold Regions in China from 1961 to 2019**

**Yumeng Wang 1, Jingyan Ma 2,\*, Lijuan Zhang 1, Yutao Huang 1, Xihui Guo 1, Yiping Yang 1, Enbo Zhao 1, Yufeng Zhao 1, Yue Chu 1, Meiyi Jiang 1 and Nan Wang 1**

**Abstract:** In this study, on the basis of the temperature data collected at 612 meteorological stations in China from 1961 to 2019, cold regions were defined using three indicators: an average temperature of <−3.0 ◦C during the coldest month; less than five months with an average temperature of >10 ◦C; and an annual average temperature of ≤5 ◦C. Spatial interpolation, spatial superposition, a trend analysis, and a spatial similarity analysis were used to obtain the spatial distribution of the cold regions in China from 1961 to 2019. Then, the areas of the cold regions and the spatial change characteristics were analyzed. The results reveal that the average area of the cold regions in China from 1961 to 2019 was about 427.70 × 104 km2, accounting for about 44.5% of the total land area. The rate of change of the area of the cold regions from 1961 to 2019 was −14.272 × 104 km2/10 a, exhibiting a very significant decreasing trend. On the basis of the changes between 1991–2019 and 1961–1990, the area of China's cold regions decreased by 49.32 × 104 km2. The findings of this study provide references for studying changes in the natural environment due to climate change, as well as for studying changes on a global scale.

**Keywords:** cold region of China; spatiotemporal distribution; spatiotemporal variation; 1961–2019

#### **1. Introduction**

Cold regions are a valuable part of the earth system [1], and they are generally defined as areas with low temperatures [2] and with the presence of ice and snow for at least part of the year. Therefore, the ice and snow resources in cold regions are very rich and are very important to human production activities [3]. They are an essential and lasting factor that affects the natural and social systems on the Earth's surface [4]. Cold regions are very sensitive to climate change [5], and climate change has changed the distribution of the cold regions around the world [6–10]. The sixth report of the IPCC pointed out that the global annual average surface temperature has risen by 1.09 ◦C in the past 100 years, that the climate warming trend in China was much higher than the global average [11,12], and that the area and spatiotemporal distribution characteristics of China's cold regions are undergoing significant changes. However, at present, little attention has been paid to the temporal and spatial changes in the cold regions in China. Gaining an accurate understanding of the distribution of, and variations in, cold regions has important theoretical and practical significance for engineering projects in cold regions, industrial and agricultural production, as well as for the rational development and utilization of the water, ice, and snow resources in cold regions [13,14].

Different methods have been proposed for the classification of cold regions. Koppen et al. [15] was the first to propose a division index for cold regions. Two indices were used

**Citation:** Wang, Y.; Ma, J.; Zhang, L.; Huang, Y.; Guo, X.; Yang, Y.; Zhao, E.; Zhao, Y.; Chu, Y.; Jiang, M.; et al. Spatial Distribution of, and Variations in, Cold Regions in China from 1961 to 2019. *Sustainability* **2022**, *14*, 465. https://doi.org/10.3390/su14010465

Academic Editors: Xiaodong Yan, Jia Yang and Shaofei Jin

Received: 9 December 2021 Accepted: 29 December 2021 Published: 2 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to classify the cold regions in Canada: (1) The average temperature of the coldest month is ≤−3.0 ◦C; and (2) The number of months with an average monthly temperature of >10 ◦C is less than four. Gerdel et al. [16] suggested dividing the Canadian cold regions on the basis of the criterion of an annual average temperature of 0 ◦C and lower. However, Wilson et al. [17] reported that there were some problems with only considering the temperature factor and proposed the use of both temperature and precipitation for cold region division. Hamelin et al. [18] proposed 10 indicators to divide the Canadian cold regions. For the classification of the cold regions in China, Yang et al. [19] proposed climate indicators on the basis of the above studies. The indicators included: (1) The average temperature in the coldest month is below −3 ◦C; (2) The number of months with an average monthly temperature of above 10 ◦C is less than four; (3) The freezing period of rivers and lakes is more than 100 days, and more than 50% of the precipitation is solid precipitation; (4) The number of months with an average monthly temperature of >10 ◦C is no greater than five; (5) The average temperature in October and April is below 0 ◦C; (6) The annual average temperature does not exceed 5 ◦C; (7) The number of days with an average daily temperature of >10 ◦C is less than 150; (8) The accumulated temperature is 500–1000 ◦C; (9) The percentage of solid precipitation is greater than 30%; and (10) The average annual number of snow cover days is >30 days. They divided the cold regions in China on the basis of these 10 indicators [3]. Their results show that China's cold regions were mainly distributed in four main regions: (1) Gansu, Qinghai, and Xinjiang; (2) Tibet, Aba, and Ganzi, in western Sichuan, northern Yunnan, the Yulong Mountains, and the north part of the Gaoligong Mountains in Yunnan; (3) Northeast and northwest Heilongjiang; and (4) The northeastern part of Inner Mongolia, except for the desert areas in the Junggar Basin, the Tarim Basin, and the northern part of Heihe [3,19]. The data reveals that the cold regions accounted for about 43% of the land area of China. To consolidate the various indicators proposed by Yang et al. [19], Chen et al. [20] proposed three indicators: the average temperature in the coldest month is <−3.0 ◦C; the number of months with an average temperature of >10 ◦C is less than five; and the annual average temperature is ≤5 ◦C. Then, using these three indicators, Chen et al. [20] created a spatial distribution map of the cold regions in China from 1961 to 1998 based on observation data collected four times a day at 571 stations in China from 1961 to 1998. They reported that China's cold regions were mainly distributed in the Greater Khingan Mountains, in the Changbai Mountains, on the Sanjiang Plain in northeastern China, in the Hexi Corridor, in most of the mountainous areas in Xinjiang, and on the Qinghai-Tibetan Plateau. The cold regions covered an area of 417.4 × 104 km2, which accounted for about 43.5% of the land area. Their results provide an important basis for studying China's permafrost, glaciers, stable seasonal snow, climatic divisions, and vegetation divisions, and these results are used to this day. However, 23 years have passed since Chen's study period ended (i.e., 1998). Since the beginning of the 21st century, the average annual temperature in China has risen by about 0.63 ◦C [21]. Therefore, it is time to revisit China's cold region divisions and to further analyze the spatial evolution characteristics of the cold regions within the context of global warming [22–24]. Currently, such studies are lacking.

In the present study, on the basis of the three indicators proposed by Chen et al. [20] and observation data collected at 612 meteorological stations in China from 1961–2019, the cold regions in China were identified, and the spatial distribution of the cold regions was analyzed in order to provide scientific references for development in the cold regions in China, and for studying the impacts of climate change in these cold regions.

#### **2. Data and Methods**

September 2020). All the 612 stations have complete data availability and there are no gaps. The distribution of the meteorological stations is shown in Figure 1.

**Figure 1.** Distribution of meteorological stations in China.

2. Elevation data: The digital elevation model (DEM) data were obtained from the cloud geospatial database (http://www.gscloud.cn/ accessed on 20 September 2020), with a resolution of 1 km × 1 km.

#### *2.2. Methods*

#### 2.2.1. Spatial Interpolation of Meteorological Data

The distribution of the meteorological stations in China is uneven. For example, in remote areas such as the Qinghai-Tibetan Plateau, there are only a limited number of stations. In order to obtain high-resolution meteorological data, spatial interpolation is often used. There are many available methods of interpolating meteorological data. Of these methods, the Kriging method is commonly used. The Kriging method performs interpolation using the distribution of the meteorological elements of nearby stations, on the basis of the principle of covariance. However, since the number of meteorological stations in cold regions is small and the spatial resolution is very low, the Kriging method is not suitable. Liu et al. [25] proposed an elevation-based spatial interpolation method, i.e., ANUSPLIN. In addition to ordinary spline independent variables, this method introduces linear covariate submodels, such as the relationship between the temperature and the elevation, precipitation, and coastline. The basic principle is to allow the introduction of multivariate and covariate linear submodels. The coefficients of the models can be automatically determined on the basis of the data. This method is able to smoothly process splines of more than two dimensions, and multiple influencing factors are introduced as covariates to perform the spatial interpolation of the meteorological elements. In addition, the temperature is affected by the altitude. Because only considering the relationship between the latitude and longitude and temperature compromises the accuracy, the terrain factor was introduced in this study as a covariate. The partial thin-plate smooth spline function of the elevation linear submodel was used for the spatial temperature interpolation. The calculation formula is as follows:

$$Z\_i = f(\mathbf{x}\_i) + b^T y\_i + e\_i(i = 1, 2, \dots, n) \tag{1}$$

where *Zi* is the dependent variable at point *i*; *xi* is a d-dimensional spline independent variable; *f* is the unknown smoothing function to be estimated; *yi* is a p-dimensional independent covariate; *b* is the p-dimensional coefficient of *yi*; *ei* is the independent variable random error with an expected value of 0 and a variance of *wiσ2*; and *wi* is the known local relative coefficient of variation, of which *σ2* is taken as the weight, and is the error variance, which is an unknown constant for all data points. When interpolating temperature, this

paper uses a three-variable local thin-disk smooth spline function (longitude and latitude are independent variables, and altitude is a covariate), and the number of splines is set to 2.

#### 2.2.2. Trend Analysis Method

In order to reflect the change trend of the cold area from 1961 to 2019, linear trend rate estimation was adopted:

$$\mathbf{y} = \mathbf{a}\mathbf{x} + \mathbf{b} \tag{2}$$

In the formula, y is the area of the cold area; x is the year; a is the linear regression coefficient, reflecting the change trend of the area of the cold area, and a > 0 indicates that the area of the cold area is increasing, anda<0 indicates that the area is decreasing; b is the intercept; and a × 10 is called the "climate tendency rate", and the unit is 104 km2/10 a.

#### 2.2.3. Spatial Similarity Analysis

The kappa coefficient is generally used to determine the degree of agreement or accuracy between two images, and its calculation formula is:

$$K = \frac{P\_0 - P\_\varepsilon}{1 - P\_\varepsilon} \tag{3}$$

Among them, *P*0 is the sum of the number of samples correctly classified in each category divided by the total number of samples, which is the overall classification accuracy. Assume that the numbers of real samples in each category are *a*1, *a*2, ..., *a*c, and that the predicted numbers of samples in each category are *b*1, *b*2, ..., *b*c, and that the total number of samples is *n*. Then:

$$P\_{\varepsilon} = \frac{a\_1 \times b\_1 + a\_2 \times b\_2 + \dots + a\_c \times b\_c}{n \times n} \tag{4}$$

The kappa coefficient calculation results are −1–1, but usually the kappa falls between 0 and 1, which can be divided into five groups to indicate the different levels of consistency: 0.0–0.20, very low consistency; 0.21–0.40, general consistency; 0.41–0.60, moderate consistency; 0.61–0.80, high consistency; and 0.81~1 are almost identical [26].

#### 2.2.4. Mann–Kendall Mutation Test

In addition to the trend analysis, the MK method can also be used to test for mutation. This method is very effective for verifying a change of state from a relatively stable state to another state. For a time series, x, with *n* sample sizes, construct an order column:

$$S\_k = \sum\_{i=1}^k r\_i \ (k = 2, 3, \dots, n) \tag{5}$$

where

$$r\_i = \begin{cases} \mathbf{1}, \mathbf{x}\_i > \mathbf{x}\_j \\ \mathbf{0}, \mathbf{x}\_i \le \mathbf{x}\_j \end{cases} \qquad (j = 1, 2, \dots, i) \tag{6}$$

It can be seen that the rank sequence, *Sk*, is the cumulative number of times the value of *i* at the moment, *i*, is greater than the number of values at time, *j*. Under the assumption of the random independence of the time series, define statistics:

$$\text{LIF}\_k = \frac{\left[\mathcal{S}\_k - E(\mathcal{S}\_k)\right]}{\sqrt{Var(\mathcal{S}\_k)}} \quad (k = 1, 2, \dots, n) \tag{7}$$

where *UF*1 = 0; and *E*(*Sk*) and *Var* (*Sk*) are the mean and variance of the cumulative number, *Sk*, respectively. This value is calculated when *x*1*, x*2*, ..., xn* are independent, and when they have the same continuous distribution as:

$$E(S\_k) = \frac{n(n-1)}{4} \tag{8}$$

$$Var(S\_k) = \frac{n(n-1)(2n+5)}{72} \tag{9}$$

*UFi* is a standard normal distribution, which is a sequence calculated according to a time series (*x*) order (*x*1*, x*2,..., *xn*). Given a significance level, a, in comparison with the data in the known normal distribution table, and if *UFi* > *Ua*, then significant changes exist in the trend. This method can also be applied to the inverse sequence of the time series, and the above procedure can be repeated by *xn, xn*−1*, ..., x*1, thus making *UFk =* −*UBk,k=n*, *n* − 1,..., and *UB =* 0. Given the significance level, *α*, the two curves of *UFk* and *UBk* and the significant horizontal line are plotted on the same graph. If the values of *UFk* and *UBk* are greater than 0, then the sequence shows an upward trend, and values below 0 indicate a downward trend. When the value exceeds the critical line, this indicates that the rising or falling trend is significant. The range beyond the critical line is defined as the time zone of mutation. If the *UFk* and *UBk* curves appear on an intersection point, and the intersection point is between the critical line, then the intersection point corresponds to the time the mutation begins. More detailed descriptions of this method are introduced in [27].

#### **3. Results and Analysis**

#### *3.1. Time Series Changes in the Area of Cold Regions in China from 1961 to 2019*

Once the data layers corresponding to the three indicators were obtained, the three layers were superimposed to extract the overlapping area, which is the spatial distribution map of the cold regions in China from 1961 to 2019. The area of the cold regions from 1961 to 2019 was calculated, and its variation with time was obtained (Figure 2a). It can be seen that the average area of the cold regions in China from 1961 to 2019 was 427.70 × 104 km2. The largest area of the cold regions occurred in 1969 (485.92 × 104 km2), and the smallest area occurred in 2007 (368.80 × 104 km2). The coefficient of variation was 0.07, indicating that the interannual variation was about 7%.

**Figure 2.** (**a**) Interannual changes; (**b**) MK test; and (**c**) Anomaly in the area of cold regions in China from 1961 to 2019.

The rate of change of the area of the cold regions in China from 1961 to 2019 was −14.272 × 104 km2/10 a, suggesting a significant decreasing trend (*p* < 0.01). Between 1961 and 2019, the area of the cold regions decreased by 84.20 × 104 km2. The Mann–Kendall (MK) test results (Figure 2b) reveal that there have been no significant sudden changes in the area of the cold regions since 1961. However, on the basis of the area anomaly map for 1961 to 2019 (Figure 2c), the area of the cold regions in China experienced a significant turning point around 1987. From 1961 to 1987, the average area of the cold regions was 453.64 × 104 km2, and from 1988 to 2019, the area decreased to 405.81 × 104 km2, i.e., a reduction of 47.83 × 104 km2. The results of the analysis of variance revealed that there was a significant difference in the area of the cold regions in China between 1961–1987 and 1988–2019 (*p* < 0.05). Thus, the area of the cold regions in China entered a declining stage in 1987. Furthermore, the rate of change of the area of the cold regions from 1961 to 1987 was −5.697 × 104 km2/10 a, indicating an insignificant decreasing trend; whereas the rate of change from 1988 to 2019 increased to −11.688 × 104 km2/10 a, and this change was significant (*p* < 0.01). In conclusion, the area of the cold regions in China has not only entered a relatively low value period since 1987, but it has also decreased significantly.

#### *3.2. Spatial Distribution of the Cold Regions in China from 1961 to 2019*

Figure 3 shows the spatial distribution of the cold regions in China from 1961 to 2019, including northeastern China, northern China, northwestern China, and southwestern China. In order to analyze the spatial distribution of China's cold regions, the cold regions were analyzed according to the administrative divisions. The cold regions were distributed in 14 provinces and autonomous regions, including Heilongjiang, Jilin, Liaoning, the Inner Mongolia Autonomous Region, Hebei, Shanxi, Shaanxi, Gansu, Ningxia, Qinghai, Sichuan, Yunnan provinces, the Xinjiang Autonomous Region, and the Tibet Autonomous Region. If you superimpose the boundaries of each administrative region on the cold area spatial distribution map calculated in this paper, using the ArcGIS Tabulate Area function, you can get the cold area of each administrative area. The areas of the cold regions in the provinces and autonomous regions are listed in Table 1, from high to low. It can be seen that the Tibet Autonomous Region contained the largest area of cold regions (105.06 × 104 km2), which was significantly larger than that of the other provinces and autonomous regions. Moreover, the areas of the cold regions in the Inner Mongolia Autonomous Region, the Xinjiang Autonomous Region, Qinghai, and Heilongjiang were also large, ranging from 54.09 × 104 km2 to 75.04 × 104 km2, which were significantly larger than those in Sichuan, Gansu, Jilin, and other provinces. By comparison, the cold regions in Hebei, Yunnan, Liaoning, the Ningxia Hui Autonomous Region, and Shaanxi were relatively small. Shaanxi had the smallest area of 0.32 × 104 km2. The entire Tibet Autonomous Region was a cold region, except for the southern part of the Nyainqen Tanglha Mountains. In the Inner Mongolia Autonomous Region, the cold regions were mainly distributed in the northeastern part of the Inner Mongolia Plateau, and in the Yinshan, Langshan, Daqingshan, and Helan mountains. In the Xinjiang Autonomous Region, the cold regions were mainly located in the Bogda Mountains in the northern part of Hami, in the Tianshan Mountains, in the Bolokonu Mountains west of Urumqi, in the Harke Mountains in the southern part of the Tianshan Mountains, and in the Altai Mountains on the Sino-Mongolian border. All of the areas in Qinghai, except for the Qaidam Basin, were cold regions. The cold regions in Heilongjiang included the Greater Khingan Mountains, the Lesser Khingan Mountains, the Changbai Mountains, the Sanjiang Plain, and the Songnen Plain. The cold regions in Sichuan Province included the high mountains in northwestern Sichuan, such as the Daxue Mountains and the Qionglai Mountains. The cold regions in Gansu Province were mainly distributed in the Qilian Mountains in the west, in the Minshan Mountains in the Longnan area, in the northern mountains in the Hexi Corridor, in the Longzhong Plateau area, and in the Mazong Mountains in the northern part of Gansu Province. The cold regions in Jilin Province were mainly located in the Changbai Mountains in the eastern part of Jilin Province and the areas east of the Hada Mountains. The cold regions in Hebei were mainly distributed in Hengshan and Yanshan in the northwest. The cold regions in Shanxi Province were mainly distributed in the high mountains, including in the Taihang Mountains, the Lvliang Mountains, the Wutai Mountains, and the Heng Mountains. The cold regions in Yunnan Province were scattered in the southern part of the Hengduan Mountains in the northwest. The cold regions in the Ningxia Hui Autonomous Region were mainly distributed in the northern part of the Liupan Mountains. In Shaanxi, the cold regions were scattered in the middle of the Qinling Mountains.

**Figure 3.** Spatial distribution of cold regions in China from 1961 to 2019.

**Table 1.** Area of cold regions in each administrative division from 1961 to 2019 (×104 km2).

#### *3.3. Spatial Variations in the Cold Regions in China from 1961 to 2019*

The World Meteorological Organization (WMO) defines the climate reference value (i.e., climatic state) as the average value of a certain meteorological element for 30 years, which is considered to be able to represent the climate of a location [28]. Thus, the 30-year average value has been used as the reference climate state in scientific research on climate and climate change. In this study, the periods of 1961–1990 and 1991–2019 were defined as two climatic states. The statistics show that the areas of the cold regions in China in 1961–1990 and 1991–2019 were 453.18 × 104 km2 and 403.86 × 104 km2, respectively. It can be seen that, against the background of global warming, the area of China's cold regions has been decreasing (i.e., by 49.32 × 104 km2). Table 2 shows the reduction in the areas of the cold regions in each administrative division. The results show that the areas of the cold regions in all 14 administrative regions decreased. Specifically, the reductions of 12.23 × 104 km2, 10.62 × 104 km2, and 8.09 × 104 km2 in Inner Mongolia, Xinjiang, and Jilin Province, respectively, were relatively large. The areas of the cold regions in Yunnan, the Ningxia Hui Autonomous Region, and Shaanxi Province decreased slightly, by 0.44 × 104 km2, 0.38 × 104 km2, and 0.27 × 104 km2, respectively.

**Table 2.** Reduction in areas of cold regions in each administrative division (×104 km2).

On the basis of the spatial distribution maps of the cold regions in 1961–1990 (Figure 4a) and 1991–2019 (Figure 4b), the kappa values of the spatial distribution maps of the two climatic states were both 0.934, indicating excellent consistency. That is, the spatial distribution of China's cold regions did not change significantly, but there were some differences. Figure 4c shows the difference in the spatial distributions of the cold regions in China in 1991–2019 and 1961–1990. From Figure 4c, the cold regions that disappeared in the second 30-year period (1991–2019) were mainly concentrated in the central and northern parts of Jilin Province, in the central-northern parts of the Inner Mongolia Autonomous Region, in the Altai region in Xinjiang, and in the central part of the Qaidam Basin. Scattered reductions also occurred in the Hengshan, Wutai, and Yunzhong mountains in the northern part of Shanxi Province, in the southern part of the Qinling Mountains in Shaanxi Province, in the southeastern part of the Tibet Autonomous Region, in Beishan, on the Longzhong Plateau in Gansu Province, in the Qionglai Mountains and Jiajin Mountains in Sichuan Province, in the Longgang Mountains and Qianshan Mountains in Liaoning Province, in the Taihang Mountains in Hebei Province, in Southern Duanyunling in Yunnan Province, and in the northern Liupan Mountains in the Ningxia Hui Autonomous Region.

**Figure 4.** Spatial distributions of cold regions: (**a**) 1961–1990; (**b**) 1991–2019; and (**c**) Difference.

#### **4. Discussion**

Currently, the spatial distribution of the cold regions in China is primarily based on the results reported by Chen et al. [20]. However, they delineated the spatial distribution of China's cold regions from 1961 to 1998 without considering the variation characteristics of the cold regions. In this study, we used the observation data collected at meteorological stations in China from 1961 to 2019 to analyze the spatial distribution of China's cold regions, and we analyzed the time series changes in the area of the cold regions and the variations in their spatial distribution. The analysis of the impact of climate change on cold regions is of great scientific significance and practical value for development in the cold regions in China, and for the rational development and utilization of the water, ice, and snow resources in these cold regions.

The data period of Chen et al. [20] was from 1961 to 1998, and the data used was from 571 meteorological stations in China. In this study, the data period was from 1961 to 2019, and the data used was from 612 meteorological stations in China. In order to compare our results with those of Chen et al. [20], the cold region distribution map of Chen et al. [20] (Figure 5b) was vectorized, and the area of the distribution map was compared with that of this study (Figure 5a). The comparison results are shown in Figure 5c. The cold area obtained in this study was 427.70 × 104 km2, accounting for 44.5% of the total area of China. The cold area obtained by Chen et al. [20] was 417 × 104 km2, accounting for 43.5% of the total area. Thus, there was a difference of 10.3 × 104 km2. As can be seen from Figure 5c, the results of this study reveal that the areas of the cold regions increased in the Qianshan and Longgang mountains in the northeastern part of Liaoning Province, in the Hengshan, Wutai, and Yunzhong mountains in Shanxi Province, in the Altay region in Xinjiang, in the Tianshan, Altun, Qilian, and Daxue Mountains, and in the Duanyun Mountains and the Mangkang Mountains on the southern Qinghai-Tibetan Plateau. In contrast, the areas of the cold regions decreased in the southern Greater Khingan Mountains, on the southern Songnen Plain, in Shanding Hural in the Inner Mongolia Autonomous Region, and in

the Qaidam Basin in Qinghai. Thus, the cold regions identified in this study are different from those identified by Chen et al. [20] because these studies were based on different research periods.

**Figure 5.** Comparison of spatial analysis results of cold regions in China: (**a**) 1961–2019; (**b**) The cold areas obtained by Chen et al.; and (**c**) Difference.

In order to compare our results with the research results of Chen Rensheng, this paper recalculated the distribution of cold regions in China from 1961 to 1998, as shown in the Figure 6 below. It can be seen that the spatial distribution map of China's cold regions calculated in this paper has higher spatial resolution than the research results of Chen Rensheng, so the statistical area should be more accurate. In addition, different spatial interpolation methods were used in the two studies, which will also make the statistical cold area different from the research results of Chen Rensheng. The main difference between the two results is due to the different time scales. More importantly, the research period of this paper is 1961–2019, which is 23 years longer than that of 1961–1998. With the extension of the research period, the spatial superposition area of the three indicators for dividing the cold region will also increase. Therefore, the area of China's cold region will increase from 1961 to 2019.

**Figure 6.** Spatial distribution of cold regions in China from 1961 to 1998.

On the basis of the area anomaly map and the results of the 5-year moving average, the year when the area anomaly became negative was 1987. Since 1987, the area of the cold regions has decreased significantly. According to the fifth assessment report of the International Panel on Climate Change (IPCC), the temperature increase became more significant after the 1980s [29]. Moreover, on the basis of the temperature data collected at 349 meteorological stations in China from 1953 to 2012, Liu et al. [30] found that the temperature increase in China during this period was significant, and that most of the years with abrupt changes were after 1986. Before the mid-1980s, the temperature in China fluctuated within a relatively small range. However, since then, the temperature has exhibited a significant upward trend. On the basis of the above discussion, it can be concluded that China's surface temperature underwent major changes in the mid-1980s. Therefore, it is concluded that the temperature change in the mid-1980s was the main influencing factor of the change in the area of the cold regions in China.

By comparing the areas of the cold regions in China under two climatic states, it was found that the area of the cold regions in northeastern China decreased by 13.62 × 104 km2, and that this was the region with the largest decrease. The results of many studies have shown that northeastern China is sensitive to climate change, and that it has experienced a significant temperature increase. For example, Chen et al. [31] studied climate change in China from 1951 to 1995 and they report that the temperature increases in China mainly occurred north of 35◦ N, with the largest temperature increase occurring in northern Heilongjiang. Haiying et al. [32] analyzed the spatial and temporal characteristics of the climate change in China from 1900 to 2000, and they found that the overall climate had been cooling since 1950, but that warming occurred in the northeastern, northern, and northwestern regions. The annual average temperature in northeastern China increased by about 1 ◦C during 1900–2000, especially from 1981 to 1998, when there was a steep jump. In addition, Liang et al. [33] analyzed the temperature characteristics in northern China from 1951 to 2014. They also found that the northeastern region had the largest rate of temperature increase, and this was also the region with the fastest temperature increase after the sudden change in temperature across the country. Therefore, the northeastern region had the largest decrease in the area of the cold regions.

The range of China's cold regions defined in this article is mainly distributed in northern Xinjiang, northeastern China, and on the Qinghai-Tibet Plateau, and this is the main distribution area of glaciers, frozen soil, and snow [34]. According to the results of this article, the area of China's cold regions is decreasing, and studies have shown that, with climate warming, glaciers shrink, the permafrost, as a whole, degenerates, snow cover shrinks, glaciers melt and shrink faster, and meltwater increases year by year. The changes in the temporal and spatial distributions of the water resources and water cycle processes caused by changes in glaciers will undoubtedly have a profound impact on the social and economic development of cold regions [35]. For example, they could lead to an increase in the glacier meltwater runoff, causing thermal melt slump, thermal melt subsidence, and other permafrost thermal melt disasters. As glacier retreat intensifies, the amount of meltwater increases, and glacier floods and glacial debris flow disasters increase with the increase in the glacier meltwater runoff [36]. Therefore, the reduction in the area of the cold region has a great impact on the ice layer, on ecosystems, and on human activities.

#### **5. Conclusions**

From 1961 to 2019, the area of the cold regions in China was about 427.70 × 104 km2, accounting for about 44.5% of China's total land area. The rate of change in the area of the cold regions was −14.272 × 104 km2/10 a, exhibiting a very significant decreasing trend. The area of the cold regions of China decreased significantly after 1987.

The cold regions in China were mainly distributed in the Greater Khingan Mountains, in the Changbai Mountains, on the Sanjiang Plain in northeastern China, in the central part of the Inner Mongolia Plateau, and in most of the mountainous areas in Xinjiang, and on the Qinghai-Tibetan Plateau. Among the 14 provinces and autonomous regions in China, the area of the cold regions in the Tibet Autonomous Region was the largest (105.06 × 104 km2), and the area of the cold regions in Shaanxi Province was the smallest (0.32 × 104 km2).

The areas of the cold regions in China during 1961–1990 and 1991–2019 were 453.18 × 104 km2 and 403.86 × 104 km2, respectively. The difference between the two climatic states was 49.32 × 104 km2. The largest decrease in the area of the cold regions occurred in the Inner Mongolia Autonomous Region (12.23 × 104 km2). The difference in the areas of the cold regions during the two periods was mainly distributed in northeastern China, in Xinjiang, and in Inner Mongolia.

**Author Contributions:** Y.W. analyzed the data and drafted the manuscript; L.Z., J.M., and Y.H. completed the manuscript and made major revisions; X.G. and Y.Y. completed the data visualization; E.Z. and Y.Z. downloaded data and searched references; Y.C., N.W., and M.J. checked and proofread the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant No. 41771067, and supported by the Key Project of Natural Science Foundation of Heilongjiang Province, No. ZD2020D002.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

#### *Article* **Snow-Disaster Risk Zoning and Assessment in Heilongjiang Province**

**Hao Li, Wenshuang Xi, Lijuan Zhang \* and Shuying Zang \***

Heilongjiang Province Key Laboratory of Geographical Environment Monitoring and Spatial Information Service in Cold Regions, Harbin Normal University, Harbin 150025, China<EMAIL_ADDRESS>(H.L<EMAIL_ADDRESS>(W.X.)

**\*** Correspondence<EMAIL_ADDRESS>(L.Z<EMAIL_ADDRESS>(S.Z.)

**Abstract:** Heilongjiang Province is located in the northeast region of China, with the country's highest latitude. It has long and cold winters, and a temperate monsoon climate. Its unique geographic location and climatic conditions make it the second largest stable snow-covered region in China. The winter snow period starts in October and ends in April of the following year. Therefore, the long-term accumulation of snow causes road obstructions and low-temperature frost damage, which seriously affects local economic development and human safety. This study adopts snow parameters (e.g., snow depth and snow-cover period), natural environmental factors (e.g., elevation and slope), and socioeconomic factors (e.g., gross domestic product and light index). On the basis of the disaster risk assessment theory, we constructed a disaster risk index from four aspects (i.e., disaster risk, susceptibility, vulnerability, and disaster prevention and mitigation). Then, we performed snowdisaster risk zoning and an assessment in Heilongjiang Province. The main findings are as follows: the snow-disaster risk in the northern and eastern regions of Heilongjiang Province was high; the central and northern regions were highly sensitive to disasters; the main urban areas were highly vulnerable; and the economically developed regions had strong disaster prevention and mitigation capabilities. Overall, the spatial distribution of the snow-disaster risk followed a decreasing trend from east to west. High-risk areas were distributed in the east and northwest (covering 34.3% of the entire Heilongjiang Province area); medium-risk areas were distributed in the north and center (accounting for 45.2% of the entire Heilongjiang Province area); and low-risk areas were concentrated in the west (constituting 20.5% of the entire Heilongjiang Province area).

**Keywords:** snow disaster; risk assessment; risk zoning; Heilongjiang Province

#### **1. Introduction**

Since the 1990s, the focus of emergency management has gradually turned from emergency rescue and postevent recovery and reconstruction, to prepreventive preparations [1]. The United States, Britain, Germany, France, Japan, and other countries have promoted major disaster risk assessments. In 2004, the United Nations International Strategy for Disaster Reduction (UNISDR) and the United Nations Development Programme (UNDP) released two reports, respectively titled: "Living with Risk: A Global Review of Disaster Reduction Initiatives", and "Disaster risk reduction: a development concern". In March 2015, the "Sendai Framework for Disaster Risk Reduction 2015–2030", finally adopted by the Third World Conference on Disaster Risk Reduction, pointed out that two of the four priorities for disaster reduction are "Understanding disaster risk" and "Strengthening disaster risk governance to manage disaster risk" [2,3]. A snow disaster is large-scale snow, caused by heavy snowfall, which seriously affects the survival and health of humans and livestock. It is a meteorological disaster that can affect and damage traffic, communications, agriculture, and electricity [4]. In 1977, the snow disaster in Xilingol (Inner Mongolia) caused more than 70% of livestock deaths [5]. From 2000 to 2012, the agriculture sector of Liaoning Province suffered 20 large-scale disasters resulting from snowstorms [6]. In 2008,

**Citation:** Li, H.; Xi, W.; Zhang, L.; Zang, S. Snow-Disaster Risk Zoning and Assessment in Heilongjiang Province. *Sustainability* **2021**, *13*, 14010. https://doi.org/10.3390/ su132414010

Academic Editors: Jia Yang, Xiaodong Yan and Shaofei Jin

Received: 28 October 2021 Accepted: 14 December 2021 Published: 18 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the 100-year snow disaster swept over half of China, and affected circuits, communications, the water supply, and heating, to different degrees. According to statistics, the snow disaster caused 129 deaths; the emergency resettlement of 1.66 million people; 485,000 house collapses; and 178 million mu of crop damage; and resulted in direct economic losses of CNY 151.65 billion. The heavy snow in Heilongjiang Province in 2007 caused 754 houses to collapse. Local snowstorms affected more than 900 people in Suifenhe City and other places, with direct economic losses of nearly CNY 100 million [7]. It is evident that snow disasters have had a huge impact on human society. Therefore, snow-disaster risk assessment provides an important theoretical basis for predisaster preparations, as well as for scientifically setting the meteorological disaster prevention standards for various regional infrastructures, which is more conducive to promoting the construction of resilient cities and villages, and to fundamentally improving the ability to resist snow disasters.

Snow disasters occur in various forms, such as snow melting leading to avalanches and ice floods, less (more) snow leading to black (white) disasters, and abnormal snowfall in windy weather leading to snowstorms. The existing studies on snow-disaster risk mainly focus on three aspects: (1) Studies on single snow-disaster risks. Some scholars have evaluated the avalanche risk. Schmitt et al. assessed the avalanche risk on the Alps on the basis of topographical parameters, such as slope, aspect and elevation, forestry-related variables, and rocks [8]. Seliverstov et al. [9] carried out avalanche risk zoning in Russia on the basis of the recurrence interval of avalanches (avalanche frequency), the percentage of the whole investigated territory that is occupied by avalanche-prone areas, the duration of the avalanche danger period, the probability of a person's stay in an avalanche-prone area for 1 day (24 h) and for 1 year, and the total population of the area and its density. Cappabianca et al. [10] presented an avalanche risk estimation procedure that combines a statistical analysis of the snowfall record, iterative simulations of avalanche dynamics, and empirically based vulnerability relations. Germain et al. [11] reconstructed past avalanche events in the north of the Gaspe Peninsula on the basis of tree rings and assessed the avalanche risk. In 2016, Germain also analyzed the avalanche risk from natural factors, the population, and environmental factors for avalanche disaster in northern Quebec, eastern Canada [12]. Some scholars have studied the risk of snowstorm disasters. Zhang et al. [13] calculated the probability of each level of the snowstorm information diffusion theory on the basis of the daily snow precipitation in 63 cities and counties in Heilongjiang Province, China, from 1961 to 2015, and established a hazard index model by using the snowstorm probability and the amount of snow precipitation. Next, a hazard assessment and the regionalization of snowstorms was performed for Heilongjiang Province from 1961 to 2015, and it was proposed that the high-risk areas increased by 30.7% from the 1960s to the 2010s, as opposed to the 38.9% reduction in the low-risk areas. Liu et al. [14] constructed the snowstorm disaster risk index on the basis of the frequency of snowstorms and assessed and regionalized the snowstorm risk in northeast China in the future (2020–2099). They showed that, under the RCP4.5 and RCP8.5 scenarios, the areas of the high-risk and lightrisk regions would increase, while the areas of the low- and medium-risk regions would decrease. (2) Studies on the impact of snow disasters. Snow disasters have a major impact on agriculture, animal husbandry, and road transportation. Barbolini et al. [15], on the basis of avalanche accidents that occurred during outdoor winter activities over the Italian Alps in the last 20 years, proposed a vulnerability relation for people directly exposed to avalanches. Sinickas et al. [16] discussed the effects of the long-term changes in the avalanche occurrence rates in terms of consequences and vulnerability. Casteller et al. [17] report on the relationship between Nothofagus broadleaf forests and avalanche runout distances. Gao [18] constructed an agricultural risk estimation index and divided the risks of China's ice and snow on agriculture. Their results propose that the regions with high snow and freezing occurrences are located in northwestern China, and that the regions with the high-loss areas are located in the coastal and southeastern parts. Sa et al. [19] established a snow-cover index using passive microwave remote sensing data from 1978 to 2012, and evaluated the risk factors of disaster in the pastoral areas of the Inner Mongolia

Grasslands. Liu et al. [20] constructed the snow-disaster index and the vulnerability index by using the highway network, social and economic development data, meteorology data, and ArcGIS software of Guoluo Prefecture, Qinghai Province, combined with the analytic hierarchy process (AHP) and cluster analysis, to calculate the snow-road-disaster risk index of Guoluo Prefecture. Qi et al. [21] proposed the risk assessment index system of highway snow cover on the basis of a theoretical model and put forward a distribution of trunk highway snow-disaster risks in Shaanxi. Leone et al. [22] evaluated the impact of avalanches on road traffic in three Alpine departments (Alpes-de-Haute-Provence, Hautes-Alpes, and Alpes-Maritimes). Snow disasters can also amplify the effects of other disasters, such as earthquakes [23]. (3) Studies on the comprehensive risk regionalization of snow disasters. Wang et al. [24] further analyzed the formation mechanism of snow disasters (SDs), and constructed the integrated risk index (IRI) of a snow disaster from the aspects of historical disasters, snowfall events, disaster-formation environments, livestock overload, and livestock vulnerability and adaptability. The regions with high IRIs are mainly concentrated in the middle, east, and southwest of the Qinghai-Tibet Plateau, and appear as a contiguous risk belt from northeast to southwest. Liu et al. [25] conducted a comprehensive analysis of the 18 indexes of snow disaster on the Qinghai-Tibet Plateau, based on the hazard harmfulness data collected from historical records and the data collected from the entities affected by this hazard in 2010, and classified the snow-disaster hazards. Gao et al. [26], on the basis of the logistic regression of 33 snow-disaster events in Qinghai province, with the maximum snow depths, snow-cover days (SCDs), slopes, annual average temperatures, and per capita gross domestic product (GDP), constructed a potential risk assessment model to regionalize the snow-disaster risk in Qinghai Province. Park et al. [27] selected the pressure index (PI), the state index (SI), and the response index (RI) to assess the comprehensive risk of snow disaster in the metropolitan city of Ulsan. Insang [28] also assessed and analyzed the snow disaster risks of Daegu City, Ulsan City, Gyeongsangbuk Province, and Gangwon Province on the basis of the subindicators of three hazards, six exposures, four vulnerabilities, and five adaptive capacities.

In summary, scholars have conducted a lot of studies on snow-disaster risk and its impact. In addition to studies on the risk of a single snow disaster, studies on snowdisaster risk regionalization, from the perspectives of risk, sensitivity, vulnerability, and disaster prevention and mitigation, have gradually increased in recent years. However, the existing risk studies of snow disaster are mostly based on snow-cover data, and the risk indicators of snowstorms and wind-blown snow are seldom considered. Therefore, the risk of snow disaster in this paper comprehensively considers snow, snowstorms, and wind-blown snow. Compared with previous studies, the index selection has certain characteristics and is comprehensive. Heilongjiang Province is the highest-latitude area in China, belonging to the second largest snow-covered area in China. It is different from other snow-covered areas because of its large annual average snow reserves and an obvious interannual variability [29]. Studies have shown the increasing intensity of winter snowfall in Heilongjiang Province in recent years [30]. The years with abnormal snowfall have also increased, especially in the first decade of the 21st century. Extreme winter precipitation events in Heilongjiang Province have frequently occurred, and the precipitation values in 2002, 2003, 2009, and 2010 were over 50% more than those in normal years [31]. Heilongjiang Province is the region with the highest average annual wind speed in China [32]. The secondary disasters caused by wind and snow are greater [7], and the persistence of snow disasters is more significant than in other regions. Therefore, we selected Heilongjiang Province as the study area. According to the natural disaster risk theory model, we established a snow-disaster risk index according to the disaster risk, susceptibility, vulnerability, and the disaster prevention and mitigation capability. Combined with geographic information system (GIS) spatial analysis tools, the weighted comprehensive evaluation method, and the analytic hierarchy process, this study classifies the risk zoning of snow disasters in Heilongjiang Province. It provides a quantitative reference for Heilongjiang province from which to determine the snow-disaster prevention

standard, make a disaster prevention plan, and assess the losses after a disaster. It is also useful for carrying out snow-disaster risk assessments and regionalization-related research.

#### **2. Materials and Methodology**

#### *2.1. Study Area*

Heilongjiang Province, located in northeastern China (121◦11 –135◦05 E, 43◦26 –53◦33 N), is the most northernly territory of the country (Figure 1). It borders Russia in the north, along Heilongjiang River. Characterized by a temperate continental monsoon climate, it has long cold winters and short summers. The annual average temperature ranges from −4 ◦C to 5 ◦C, from north to south of the study area. Its total annual precipitation is around 400–650 mm, which gradually decreases from east to west. The average annual wind speeds are about 2–4 m/s. Heilongjiang Province is an important grain production base in China, with commercial grain accounting for about one-tenth of the national total.

**Figure 1.** The study area and distribution of 62 meteorological stations in Heilongjiang Province, China.

#### *2.2. Datasets*

Precipitation and wind speed data: The daily meteorological data used in this study were provided by the Heilongjiang Meteorological Data Center. The earliest meteorological station in Heilongjiang Province was built in 1951, and 83 stations were established later. They have been divided into national standard meteorological stations and general meteorological stations. According to the continuity of data and the number of stations required for analysis, we selected 62 stations with continuous and complete observation data since 1961. The spatial distribution of the stations is shown in Figure 1. The basic data used in the study include the daily temperatures and precipitation data. The time scale is from 1983 to 2020. The study period is winter, which is defined as taking place from November of the current year to February of the second year.

Elevation, slope, and other data: The elevation, slope, and lighting data were generated by visual interpretation from the Resource and Environment Science and Data Center (www.resdc.cn. Accessed on 9 February 2021), with a spatial resolution of 90 m. The terrain standard deviation was calculated on the basis of the elevation.

Normalized difference vegetation index (NDVI) data: The NDVI data were downloaded from the National Aeronautics and Space Administration (NASA) website, from November of the current year to February of the following year, with a resolution of 1 km × 1 km. The time period is from 1996 to 2020, which is the winter average of the past five years.

Nighttime light index data: This study used national stable light images. Briefly, a stable light image is an annual raster image that calibrates the average nighttime light intensity. Each image included the lights of cities, towns, and the permanent light sources of other places (excluding the influence of occasional noises, such as moonlight clouds, firelight, and oil and gas burning). The pixel DN value of the image represents the average light intensity, and its range was 0–63. The nighttime light index data were downloaded from the Geographical Information Monitoring Cloud Platform (http://www.dsac.cn/ DataProduct/Detail/201116. Accessed on 13 January 2021).

Socioeconomic factors: The socioeconomic factors (e.g., gross domestic product [GDP], per capita disposable income, and education level) used in this study were from the Heilongjiang Statistical Yearbook 1996–2020 (Accessed on 10 January 2021).

#### *2.3. Methods*

On the basis of the basic theory of disaster risk assessment, this paper constructs a theoretical model of a snow-disaster risk assessment. The trend analysis method is used to analyze the time characteristics of the risk factors; the analytic hierarchy process is used to determine the weights of the risk-influencing factors, and the comprehensive weighted evaluation method is used to form a risk index model.

#### 2.3.1. Basic Theory of Disaster Risk Assessment

According to the theory of natural disaster risk formation, meteorological disaster risk is formed by four parts: risk, susceptibility (hazard-bearing object), vulnerability (hazard-bearing object), and the disaster prevention and mitigation capability. Each factor is composed of a series of subfactors. Its expression is:

#### Disaster risk index = f (risk, susceptibility, vulnerability, disaster prevention and mitigation capability) (1)

The risk factors: All meteorological factors that may cause disasters can be called "meteorological hazards", and most of the meteorological hazards that exist in disastercausing environments are abnormalities of certain natural phenomena, temporal and spatial laws, or natural substances, or a certain abnormality in the process of energy exchange. Generally, the greater the risk of meteorological factors, the greater the risk of meteorological disasters.

Susceptibility of the disaster-bearing object: The disaster-bearing object is the target of the hazard-causing factor and is the entity that bears the disaster. The susceptibility of the hazard-bearing object is the property shown by the hazard-bearing individual exposed to the hazard-causing agent and is the result of the interaction between the hazard-causing factor and the hazard-bearing object.

Vulnerability of disaster-bearing object: Disasters can only become disasters when they act on corresponding objects, namely, humans and their social and economic activities. Specifically, a disaster-bearing object refers to all objects that may be threatened by hazards in a given dangerous area, as well as to the degree of damage or loss caused by potential hazards, which comprehensively reflects the degree of loss of meteorological disasters.

The disaster prevention and mitigation capability refers to the various management measures and countermeasures that are used to prevent and mitigate meteorological disasters, including management capabilities, disaster-reduction investment, resource preparation, etc. With the proper management measures and strong management capabilities, the smaller the potential loss that may be suffered, and the smaller the risk of meteorological disasters.

On the basis of the above theory, we built a hierarchical analysis model for snowdisaster risk assessment (Figure 2). Please refer to Sections 3.1–3.4 for the selection basis of the specific indicators.

**Figure 2.** Hierarchical model of snow-disaster risk assessment.

2.3.2. Weighted Comprehensive Evaluation Method

The weighted comprehensive evaluation method is a method for solving the "bottomup" indexes in the risk hierarchy analysis and the evaluation model [33]. This method comprehensively considers the degree of influence of each factor on the overall object and integrates the advantages and disadvantages of each specific index, using a numerical index to concentrate the data, indicating the pros and cons of the entire evaluation object. Therefore, this method is particularly suitable for the comprehensive analysis, evaluation, and optimization of technologies, strategies, or programs, and is currently one of the most commonly used calculation methods. Its expression is:

$$Y\_i = \sum\_{j=1}^{m} \lambda\_{i\bar{j}} X\_{i\bar{j}} \quad i = 1, 2, 3, 4; \quad j = 1, 2, \dots, m \tag{2}$$

In the formula, *Yi* represents the disaster risk index, and *i* represents, respesctively, the risk, susceptibility, vulnerability, and disaster prevention and mitigation capabilities; *Xij* is the factor that affects the risk, susceptibility, vulnerability, and disaster prevention and mitigation capability; and *λij* is the weight value of the risk, susceptibility, vulnerability, and disaster prevention and mitigation capability (0 ≤ *λj* ≤ 1).

The expression for the comprehensive risk index of natural disasters is:

$$\mathcal{Y} = \sum\_{i=1}^{n} W\_i \mathcal{Y}\_i \quad i = 1, 2, 3, 4 \tag{3}$$

In the formula, *Y* represents the comprehensive disaster risk index; *Yi* is the risk index, the susceptibility index, the vulnerability index, and the disaster prevention and mitigation capability index; and *Wi* is the weight value. The stronger the disaster prevention and mitigation ability, the smaller the comprehensive risk index and, thus, the "negative sign" is used.

The *λj* and *Wi* are determined by the analytic hierarchy process (see research method in Section 2.3.3 for details); each factor in the formula needs to be standardized because of different dimensions (see research method in Section 2.3.4 for details).

#### 2.3.3. Analytic Hierarchy Process

The analytic hierarchy process (AHP) is a simple method for making decisions on some more complex and vague problems, especially for those problems that are difficult to fully quantitatively analyze [34]. This paper uses the operation principle of the analytic hierarchy process and uses the 1–9 scale method, provided by Saaty, to construct the judgment matrix for the pairwise relationships of the influence factors. The pairwise comparison of all of the influence factors determines the weight of each influence factor, which avoids the result error caused by the subjectivity of the expert. The qualitative comparison scale values between the two influencing factors are shown in Table 1 below.

**Table 1.** Scale of the AHP method.

Solve the maximum eigenvector value of the judgment matrix and its corresponding eigenvector by the sum-product method and check the consistency of the matrix (the following formula). After passing, solve it by the sum-product method.

$$CI = \frac{\lambda\_{\max} - n}{n - 1} = \frac{-\sum\_{i=1}^{n} \lambda\_i}{n - 1} \tag{4}$$

$$CR = \frac{CI}{RI} < 0.1\tag{5}$$

In the formula, *CI* is the consistency index of the judgment matrix; *λmax*is the largest characteristic root of the matrix; *n* is the order of the discrimination matrix; *CR* is the random consistency index of the judgment matrix; *RI* is the average random consistency index of the discrimination matrix. The *RIs* of the risk, susceptibility, vulnerability and the disaster prevention and mitigation capability are 1.24, 1.12, 0.9, and 0.58, respectively.

According to the analytic hierarchy process, taking risk as an example, the process of determining the weight of each factor is as follows:

The hazard factors include the snowfall intensity, snow depth, snow-period length, snowfall days, blizzard volume, blizzard number, and the days with wind speeds over 8 m/s. It is believed that, according to the degree of importance, the relative hazard factors of the blizzard amount are the most important, followed by the number of blizzards, the snowfall intensity, the snow depth, the snow-cover period, the snowfall days, and the wind speeds. According to the relative importance of each factor, the judgment matrix is constructed according to Table 1. The maximum characteristic root of the matrix, *λmax*, is 6.242, *CR* = 0.038 < 0.10, and the weight result is valid. Therefore, the weights of the snow intensity, snow depth, snow period, the total blizzard volume, the total number of blizzards, and wind speeds over 8 m/s are 0.3848, 0.3054, 0.1648, 0.0718, 0.0468, and 0.0264, respectively.

In the same way, the weights of the susceptibility, vulnerability, disaster prevention and mitigation capability, and the comprehensive disaster risk factors are shown in Table 2.

**Table 2.** The weights of the risk, susceptibility, prevention and mitigation capacity, and comprehensive risk.

Note: TSD (terrain standard deviation); UNL (urban night lighting).

#### 2.3.4. Standardization

In the process of zoning, because of the different dimensions of the selected factors, the magnitude of the difference is large. For example, the length of the snow-cover period is 150 days, and the average number of snowfalls is about 10 times per year. Therefore, when calculating the hazard-factor risk index, it needs to be normalized so that the value of each factor is between 0 and 1. When assessing the hazard risk, the hazard-bearing body susceptibility, the hazard-bearing body vulnerability, the disaster prevention and mitigation capability, the relationships between the selected influencing factors and the hazard-causing factor risk, the hazard-bearing body susceptibility, the hazard-bearing body vulnerability, and the disaster prevention and mitigation capability are different. Some are the greater the number of influencing factors, the greater the hazard risk, the greater the susceptibility of the hazard-bearing body, the greater the vulnerability of the hazard-bearing body, and the greater the ability of disaster prevention and mitigation, while some factors are the opposite. Therefore, in the evaluation process, the maximum value standard or the minimum value standardization are performed for ostentation, and an example of the formula is as follows: The greater the snowfall intensity, the greater the hazard risk of the hazard factor. Therefore, the maximum value of the snowfall intensity is standardized, and Formula (6) is selected; if the highway mileage is longer, the vulnerability is smaller and, thus, the highway mileage is minimized. For the value standardization, select Formula (7).

Maximum standardization:

$$X'\_{\text{max}} = \frac{\left| X\_{ij} - X\_{\text{min}} \right|}{X\_{\text{max}} - X\_{\text{min}}} \tag{6}$$

Minimum standardization:

$$X'\_{\min} = \frac{|X\_{\max} - X\_{ij}|}{X\_{\max} - X\_{\min}} \tag{7}$$

where *Xij* is the index number of the *j*-th factor of the *x* factor; *X max* and *X min* are the dimensionalities of *Xij*; and *Xmax* and *Xmin* are the minimum and maximum values in the index sequence, respectively.

#### 2.3.5. Mann–Kendall Trend Test

In order to test the change trend of the meteorological elements, the Mann–Kendall trend test method is used to test the meteorological elements on the annual and seasonal scales [35]. The method assumes that the time data series (*x*1, *x*2, ... , *xn*) are independent, random, and uniformly distributed. For the sample sequence, the calculation equation for the test statistical variable, *S*, is:

$$S = \sum\_{i=2}^{n} \sum\_{j=1}^{i-1} \text{sign}(\mathbf{x}\_i - \mathbf{x}\_j) \tag{8}$$

when:

$$\mathbf{x}(\mathbf{x}\_{i} - \mathbf{x}\_{j}) = \left\{ \begin{array}{c} > 0 \\ = 0 \\ < 0 \end{array} \right\}, \operatorname{sign}(\mathbf{x}\_{i} - \mathbf{x}\_{j}) = \left\{ \begin{array}{c} 1 \\ 0 \\ -1 \end{array} \right\} \tag{9}$$

$$
\sigma s = \sqrt{\frac{n(n-1)(2n+5)}{18}}\tag{10}
$$

when:

$$S\left\{ \begin{array}{c} > 0\\ = 0\\ < 0 \end{array} \right\}, Z = \left\{ \begin{array}{c} \frac{S-1}{\sigma s} \\ 0 \\ \frac{S+1}{\sigma s} \end{array} \right\} \tag{11}$$

where *S* is the test statistical variable of the normal distribution; *xi* and *xj* are two series of different distributions in the same sample, where 1 ≤ *j* < *I* ≤ *n*; *σS* is the standard deviation; *n* is the total number of samples; and *Z* is the test value. If *Z* > 0, the tested time series has an upward trend, and if *Z* < 0, the tested time series has an upward and downward trend; the absolute value of *Z* is greater than or equal to 2.32, 1.64, and 1.28, which means they pass confidence. They are, respectively, 99%, 95%, and 90% significance test levels.

#### 2.3.6. The Pettitt Test

The Pettitt test is a common tool to detect a single unknown mutation point, and it is also one of the most common nonparametric test methods [36]. It can be described as giving an observation data sequence, *X* | *xt*, *t* = 1, 2, ... , n, where n is the sample size, if the sequence has a change point at τ(1 ≤ τ ≤ n − 1), defining the corresponding Pettitt statistics, *Uτ*.

$$\text{U}\,\text{I}\,\tau, n = \sum\_{i=1}^{\tau} \sum\_{j=\tau+1}^{n} \text{sgn}(\mathbf{x}\_{i} - \mathbf{x}\_{j}) = \begin{cases} -1(\mathbf{x}\_{i} - \mathbf{x}\_{j}) < 0 \\ 0(\mathbf{x}\_{i} - \mathbf{x}\_{j}) = 0 \\ +1(\mathbf{x}\_{i} - \mathbf{x}\_{j}) > 0 \end{cases} \tag{12}$$

Let *K* = max(|*Uτ*,*n*|), then the time, *T*, where *K* is located, is the possible change point position. The significance level of the corresponding change point is:

$$P = 2\exp\left[-6K^2(n^3 + n^2)\right] \tag{13}$$

If *P* ≤ 0.5, then time, *T*, is considered to be the position of the statistically significant change point.

#### **3. Results**

*3.1. Snow-Disaster Risk Analysis*

3.1.1. The Characteristics of the Risk Zoning Factors of Snow Disasters in Heilongjiang Province

In this paper, the snow depth, the length of the snow-cover period, the snowfall intensity, the total blizzard volume, the total number of blizzards, and the total days of wind speeds over 8 m/s were selected as the risk factors of snow disaster. The deeper the snow depth, the longer the snow-cover period, the stronger the snowfall intensity, the greater the blizzards, the greater the total number of blizzards, and the more days with wind speeds over 8 m/s, the greater the risk of snow disasters. Figure 3 shows the interannual variation characteristics of the snow depth, the snow-period length, the snow intensity, the total blizzard volume, the total number of blizzards, and the days with wind speeds over 8 m/s in Heilongjiang Province, from 1983 to 2020. It can be seen that the snow depth, snow intensity, total blizzard volume, and the total number of blizzards in Heilongjiang Province show a significant increase, and the M–K statistics are 2.869, 3.269, 1.968, 2.014, respectively. The length of the snow-cover period and the wind speeds showed a significant decrease, and the M–K statistical values were −3.223 and −6.43 (Table 3). It can be seen from the distribution of the box map (Figure 4) that the snow depth, the snowperiod length, the snow intensity, the total blizzard volume, the total number of blizzards,

and the days with average wind speeds over 8 m/s are 20.82 cm, 165.4 d, 1.24 cm/time, 142.0 mm, 9.9 times, and 92 d, respectively. The interannual changes in the wind speed, snow intensity, the number of blizzards, and the snow depth are relatively large, while the interannual changes in the snow depth and the blizzard volume are relatively small. It can be seen from the extreme value distribution characteristics that the maximum value of the snow intensity, the maximum value of the blizzard volume, and the number of blizzards are obviously larger than the average values, while the minimum value of the days with wind speeds over 8 m/s is obviously smaller than the average value. Using the Pettitt method to calculate the mutation year of each hazard factor, it was concluded that the snow depth, the snow-cover length, the snow intensity, and the days with wind speeds over 8 m/s days in Heilongjiang Province all have abrupt changes around the early 2000 s, and the catastrophes of the total snow cover and the total number of snowstorms occur in the middle and late parts of the 2000 s, respectively.

**Figure 3.** Time variation of risk factors.

**Figure 4.** Boxplots of risk factors.

Figure 5 shows the spatial distribution of the snow depth, the snow-period length, the snow intensity, the total blizzard volume, the total number of blizzards, and the days with wind speeds over 8 m/s in Heilongjiang Province, from 1983 to 2020. It can be seen that the snow depths of Heilongjiang Province, from 1983 to 2015, are between 7.39 and 35.26 cm, and that the snow depths of 25–35.26 cm are mainly distributed in Mohe County in the northwest, northwest of Tahe County, northwest of Huzhong District, Xunke County, Jiayin County, north of Yichun City, north of Hegang City, and north of Luobei County; snow depths of 15–25 cm are mainly distributed in most cities and counties in central and southern Heilongjiang Province; snow depths of 7–15 cm are mainly distributed in the southwest of Heilongjiang Province, including Qiqihar cities and counties, Daqing, and the southwest of Suihua. The length of the snow period in Heilongjiang Province from 1983 to 2015 is 150–200.3 days. Snow periods of more than 180 days are mainly distributed in the Daxing'anling area, in the north of Heilongjiang Province; the areas of snow periods of more than 160 days are mainly distributed in central and northern Heilongjiang Province; the areas with a snow-cover period of less than 160 days are mainly

distributed in the southern part of Heilongjiang Province. The length of the snow-cover period of the province generally increases from low latitudes to high latitudes. The snow intensity in Heilongjiang Province from 1983 to 2020 was from 1.23–2.10 cm/time. The areas with higher snow intensities are mainly located in the east of Hegang, Jiamusi, and in Shuangyashan in the east; the areas with moderate snow are mainly located in the east of Heilongjiang Province and Mohe County; and the areas with low snow intensity are mainly distributed in the west of Heilongjiang Province. The total blizzards in Heilongjiang Province from 1983 to 2020 were between 0 and 761.9 mm, and the blizzards above 200 mm were mainly distributed in Jiamusi, Hegang, Shuangyashan, Qitaihe, Jixi, and Mudanjiang, in eastern Heilongjiang Province; the areas with blizzards above 100 mm were mainly distributed in the Yichun, eastern Suihua, eastern Harbin, Mudanjiang, and Daxing'anling regions, in the central part of Heilongjiang Province; the areas with blizzards below 100 mm were mainly distributed in Heihe, Qiqihar, Daqing, Suihua, and Harbin, in the western part of Heilongjiang Province. The total number of blizzards in Heilongjiang Province from 1983 to 2020 was between 5 and 20 times. The total number of blizzards above 15 were mainly located in Jiamusi, Hegang, Shuangyashan, Qitaihe, Jixi, and Mudanjiang, in eastern Heilongjiang Province; the number of blizzards were 10 in the areas mainly distributed in Qiqihar, Daqing, Suihua, and Harbin, in the west of Heilongjiang Province. In Heilongjiang Province, from 1983 to 2020, the number of days with wind speeds over 8 m/s were 70–187.7 days, and the days with wind speeds over 8 m/s over 110 days were mainly distributed in Jiamusi, Hegang, Shuangyashan, Qitaihe, Jixi, and Mudanjiang, in eastern Heilongjiang Province. The areas with wind speeds> 8 m/s within 90–130 days are mainly distributed in the Qiqihar, Daqing, Heihe, and Daxing'anling areas, in western Heilongjiang Province; the areas with wind speeds over 8 m/s below 90 days are mainly distributed in Heihe, Suihua, and Iraq, in central Heilongjiang Province.

**Figure 5.** Distribution map of risk factors of snow disaster.

#### 3.1.2. Results of the Risk Zoning of Snow Disasters in Heilongjiang Province

The spatial superpositions of the risk factors are used to obtain the disaster risk zoning results of snow disasters in Heilongjiang Province (Figure 6a). The spatial distribution of hazards is uneven, and there are high-value areas locally. The risk index is equally divided into three levels (as shown in Table 4), and the spatial distribution of the risk classification is presented in Figure 6b. The spatial distribution of the snow-disaster risk in Heilongjiang Province followed a decreasing trend from east to west. High-risk areas accounted for 46.3% of the entire Heilongjiang area, and were distributed mainly in the eastern region of Heilongjiang (Jiamusi, Yichun, Jixi, Qitaihe, Hegang, Shuangyashan, Xunke, and Mudanjiang) and the northwestern Daxing'anling area; the medium-risk areas covered 36.8% of the entire Heilongjiang area, and were distributed primarily in the central region of Heilongjiang Province (eastern Daxing'anling, Heihe, eastern Suihua, northern Mudanjiang, and eastern Harbin); and the low-risk areas accounted for 16.9% of the total area of Heilongjiang Province, and were distributed mostly in the western region of Heilongjiang Province (Qiqihar, Daqing, western Suihua, and western Harbin).

**Figure 6.** Risk map of snow disaster in Heilongjiang Province ((**a**): stretch, (**b**): classification).

#### *3.2. Susceptibility Analysis of Snow Disaster*

3.2.1. The Characteristics of the Susceptibility Zoning Factors of Snow Disasters in Heilongjiang Province

Figure 7 shows the spatial distribution of the altitude, slope, land-use type, terrain standard deviation, and the NDVI in Heilongjiang Province. It can be seen that the altitude in Heilongjiang province is between 15 and 1651 m; the slope is between 0 and 57.4◦; the terrain standard deviation is between 0 and 174.14; and the spatial distribution of the altitude, the slope, and the terrain standard deviation are roughly the same. The areas with high levels are mainly distributed in Huzhong District, Xinlin District, Yichun City, Hailin City, Muling City, and Dongning County in Heilongjiang Province. The land-use types in Heilongjiang Province are dominated by farmland and forest. The forest land is mainly distributed in the Greater Khingan Mountains, in the Lesser Khingan Mountains, and in Mudanjiang City in the south. The NDVI values corresponding to the forest area are also higher, and the highest is 0.83. Farmland is mainly distributed on the Songnen Plain, in the west of Heilongjiang Province, and on the Sanjiang Plain in the east.

**Figure 7.** Heilongjiang Province altitude, slope, land-use type, terrain standard deviation, NDVI spatial distribution maps.

#### 3.2.2. Results of the Susceptibility Zoning of Snow Disasters in Heilongjiang Province

The weight in Table 2 is used to perform the spatial superposition of the susceptibility factors to obtain the susceptibility zoning results of snow disasters in Heilongjiang Province (Figure 8a). The spatial distribution of the snow-disaster susceptibilities in Heilongjiang Province is high in the central and northwestern regions, and is low in the eastern and western regions. The susceptibility index is equally divided into three levels (as shown in Table 4), and the spatial distribution of the susceptibility classification is presented in Figure 8b. High susceptibilities accounted for 16.3% of the entire Heilongjiang area, and were distributed mainly in the northern region (southern Mohe County, Huzhong District, Xinlin District, Yichun City, Tieli City, and Mudanjiang City); medium susceptibilities covered 35.6% of the entire Heilongjiang area, and were distributed largely in the central region of Heilongjiang Province (Heihe City, Huma County, Tahe County, Jixi City, and Yichun); and low susceptibilities accounted for 48.1% of the total area of Heilongjiang Province, and were distributed mainly in the western (Qiqihar, Daqing, Suihua, and Harbin) and eastern (Jiamusi, Shuangyashan, and Jixi) regions of Heilongjiang Province. Heilongjiang Province is located in the Northeast Plain and includes the Songnen Plain and the Sanjiang Plain. Other regions with high susceptibilities were distributed in the Daxinganling, Xiaoxinganling, and Mudanjiang mountainous regions.

#### *3.3. Vulnerability Analysis of Snow Disaster*

3.3.1. The Characteristics of the Vulnerability Zoning Factors of Snow Disasters in Heilongjiang Province

Figure 9 shows the spatial distributions of urban lights, livestock numbers, crop-sown areas, and highway mileages in Heilongjiang. It can be seen that the urban night lights are between 0 and 63, and that the high values of urban night lights are scattered in prefecturelevel cities and urban areas, while urban night lights in country-level cities are relatively low. The livestock numbers are from 498 to 665,511, and the areas with high levels are mainly distributed in Longjiang County, Shuangcheng City, and Zhaodong City, in the southwest of Heilongjiang Province. The crop-sown areas are from 214 to 342,188 km2, and

the areas with high levels are mainly distributed in Songnen Plain and Sanjiang Plain. The highway mileage is between 92.88 and 25,455, and the areas with high levels are distributed in Qiqihar city, Harbin city, Suihua city, and Heihe city, in Heilongjiang Province.

**Figure 8.** Map of snow-disaster susceptibilities in Heilongjiang Province ((**a**): stretch, (**b**): classification).

**Figure 9.** Spatial distributions of urban lights, livestock numbers, crop-sown areas, and highway mileages in Heilongjiang.

3.3.2. Results of the Vulnerability Zoning of Snow Disasters in Heilongjiang Province

The weight in Table 2 is used to perform the spatial superposition of the vulnerability factors to obtain the vulnerability zoning results of snow disasters in Heilongjiang Province (Figure 10a). The spatial distribution of the snow-disaster vulnerabilities in Heilongjiang Province is high in the western regions, and low in the central and northwestern regions. The vulnerability index is equally divided into three levels, and the spatial distribution of the vulnerability classification is presented in Figure 10b. High vulnerabilities accounted for 3.0% of the entire Heilongjiang area, and were distributed mainly in the economically developed main urban areas (Qiqihar City, Daqing City, Harbin City, and Suihua City), as well as other small parts of cities and counties; medium vulnerabilities constituted about 27.9% of the entire Heilongjiang area, and were distributed chiefly in Heihe City, Nenjiang County, Nehe City, Longjiang County, Zhaodong City, Wuchang City, Fujin City, Baoqing County, and other cities and counties of Heilongjiang Province; and low vulnerabilities accounted for 69.1% of the entire area of Heilongjiang Province, and were distributed mostly in the northern and central regions of Heilongjiang Province. As shown in the figure, high vulnerabilities were distributed mainly in regions with more urban lights and larger populations, which have been significantly affected by disasters.

**Figure 10.** Vulnerability zoning map of snow disaster in Heilongjiang Province ((**a**): stretch, (**b**): classification).

#### *3.4. Analysis of Snow-Disaster Prevention and Mitigation Capability*

3.4.1. The Characteristics of the Prevention and Mitigation Capability Zoning Factors of Snow Disasters in Heilongjiang Province

Socioeconomic factors are the main indicators that affect the disaster prevention and mitigation capability. We selected the economic indicators of the Heilongjiang Province Statistical Yearbook, from 1996 to 2020 that each county has data for, as the factors of the disaster prevention and mitigation capability, and finally selected the total GDP, the per capita disposable income, and the education level. The higher the GDP, the stronger the region's economic strength, and the stronger its ability to fight disasters. The population size of each region is different, and the per capita disposable income further illustrates the level of economic development. Therefore, the per capita disposable income was selected as the disaster prevention and mitigation ability. The level of education represents the level of awareness of the impacts of, and the defenses against, disasters. Therefore, the higher the level of education, the stronger the ability to reduce disasters.

Figure 11 shows the spatial distributions of the total GDPs, per capita disposable incomes, and education levels in Heilongjiang. It can be seen that the total GDP is between 25,538 and 15,324,071, and that the high values of the total GDPs are mainly distributed in Daqing City, Qiqihar City, and Harbin City. The per capita disposable income is between 2514 and 14,850, and the areas with high levels are mainly distributed in Daqing City, Harbin City, Dongning County, and Hailin City. The education level is between 2300 and 73,770, and the areas with high levels are the same as those with high total GDPs.

**Figure 11.** Spatial distributions of total GDPs, per capita disposable incomes, and education levels in Heilongjiang Province.

3.4.2. Results of the Prevention and Mitigation Capacity Zoning of Snow Disasters in Heilongjiang Province

The weight in Table 2 is used to perform the spatial superposition of the prevention and mitigation capacity factors to obtain the prevention and mitigation capacity zoning results of snow disaster in Heilongjiang Province (Figure 12a). The spatial distribution of the snow-disaster prevention and mitigation capacities in Heilongjiang Province is uneven and is high in the southern regions, and low in the northwestern regions. The prevention and mitigation capacity index is equally divided into three levels, and the spatial distribution of the prevention and mitigation capacity classification is presented in Figure 12b. High disaster prevention and mitigation capabilities accounted for 7.4% of the entire Heilongjiang area, and were distributed mainly in Qiqihar City, Daqing City, Zhaodong City, Hulan District, Harbin City, Acheng District, Shuangcheng City, Mudanjiang City, Jiamusi City, Suifenhe City, and Dongning County; middle disaster prevention and mitigation capabilities constituted 84.0% of the entire Heilongjiang area, and were distributed in most cities and counties in Heilongjiang Province; and low disaster prevention and mitigation capabilities covered 8.6% of the entire Heilongjiang Province area, and were distributed mainly in Sunwu County, Gannan County, Tailai County, Lindian County, Kedong County, Baiquan County, Mingshui County, Qinggang County, Wangkui County, Lanxi County, Suileng County, and Jixi City. The figure shows that the more developed the economy, the stronger the disaster prevention and mitigation capability.

**Figure 12.** Map of Heilongjiang Province's snow-disaster prevention and mitigation capacities ((**a**): stretch, (**b**): classification).

#### *3.5. Comprehensive Snow-Disaster Risk Assessment and Zoning*

The risk factors, the susceptibility factors, the vulnerability factors, and the prevention and mitigation capacity factors are superimposed, according to the weights, and the comprehensive snow disaster risk index is obtained. The comprehensive risk assessment and zoning results for snow disasters in Heilongjiang Province are shown in Figure 13a. The comprehensive risk index is equally divided into three levels (Table 4), and the spatial distribution of the comprehensive risk classification is presented in Figure 13b.The spatial distribution of the comprehensive snow disaster risk in Heilongjiang Province followed a declining trend from east to west. High risks accounted for 34.3% of the entire Heilongjiang Province area (i.e., an area of 97,000 km2), and were distributed mainly in Hulin City, the eastern part of Mishan City, Shuangyashan City, Qitaihe City, Suifenhe City, Dongning County, Muling City, Ning'an City, Hailin, Mudanjiang City, Yichun City, Jiayin County, Xunke, Hegang, Tangyuan, Jiamusi, Tieli, Mohe County, Huzhong District, Xinlin District, and western Tahe County; medium vulnerabilities covered 45.2% of the Heilongjiang area (i.e., an area of 214,000 km2), largely distributed in eastern Tahe County, Huma County, Heihe City, Nenjiang County, Sunwu County, eastern Wudalianchi, Bei'an, Qing'an, Hailun, Suihua, Mulan, Tonghe, Yilan, Luobei County, Suibin County, Tongjiang City,

Fuyuan County, Fujin, Baoqing, Youyi, Jixian, Linkou, Jixi, Jidong, Shangzhi City, and Wuchang City; and low vulnerabilities accounted for 20.5% of the Heilongjiang Province area (i.e., an area of 162,000 km2), and were distributed mostly in the southwestern region of Heilongjiang Province, including Qiqihar, Daqing, western Suihua, and western Harbin.

**Figure 13.** Comprehensive risk zoning map of Heilongjiang snow disaster ((**a**): stretch, (**b**): classification).

#### **4. Discussion**

(1) This article evaluates and zones the snow-disaster risk in Heilongjiang Province from four aspects: risk, susceptibility, vulnerability, and the disaster prevention and mitigation capability. Compared to the snow-disaster risk research at home and abroad, the analysis content of this article is complete. The zoning indicators are relatively complete, which provides a research example for related research and has a certain reference effect. It should be noted that the constituent factors of risk, susceptibility, vulnerability, and the disaster prevention and mitigation capability are not necessarily fixed, and adjustments and replacements can be made without affecting the physical meaning when conducting research in different regions. Generally speaking, six index meteorological indicators, which included the snow depth, the snow-period length, snow days, total blizzard volume, the number of blizzards, and their scores were selected for risk, which are generally easier to obtain; the altitude and land-use types that constitute susceptibility are also easier to obtain. The standard deviations of the slope and terrain are calculated on the basis of altitude data. The NDVI data uses NDVI data; and the factors affecting vulnerability and the disaster reduction capabilities are mainly derived from statistical yearbooks. Since the indicators of the statistical yearbooks in different regions are different, the data in the statistical yearbooks can be replaced when conducting research. For example, the urban lighting index was used in this article because there is no spatial distribution data of populations in the statistical yearbooks. This article uses the MODIS lighting index to replace it.

(2) In the process of zoning, not only can we obtain a spatial distribution map of the different factors that affect not only the risk, susceptibility, vulnerability, and disaster prevention and mitigation capability, but also a zoning map of the risk, susceptibility, vulnerability, disaster prevention and mitigation capability, and the comprehensive disaster risk. Therefore, this method is conducive to the comprehensive analysis of disaster risk. For example, the Sanjiang Plain, in the east of Heilongjiang Province, is a high-risk area for snow disasters, compared with Figure 14, and most of the final comprehensive snowdisaster risk is in a medium-risk area. The main reason is the low susceptibility and the low vulnerability. The comparison and influencing factors can be further analyzed as follows: Although the snow depth, intensity, blizzard volume, and the number of blizzards are relatively large and strong here, it is a plain with a low altitude, small slopes, and small topographical fluctuations. This area is, especially, a primarily agricultural area, basically

consisting of arable land, and there are almost no crops in winter, so the comprehensive risk of snow disaster is low. Therefore, the results of the snow-disaster risk zoning in this paper provide theoretical support for the government and for disaster prevention and relief departments to formulate predisaster mitigation plans, and to help with postdisaster relief decision-making.

(3) Since the selection of zoning factors in this paper will be affected by human factors, on the one hand, and will be limited by statistical yearbook data, on the other, the uncertainty of the factor selection will affect the accuracy of the zoning results. Especially when zoning a large spatial area, it is difficult to obtain a unified socioeconomic index, so it will be restricted. It is recommended to use remote sensing products as much as possible, as with the NDVI and the urban lighting index in this article, which can achieve large-scale spatial risk zoning while ensuring high spatial resolution.

In addition, many meteorological disasters, such as droughts and floods, are related to precipitation data. The precipitation is affected by small-scale weather systems, as well as by the local topography, the underlying surface, and other factors. The spatial interpolation accuracy of precipitation will be affected. The precipitation data in the published precipitation reanalysis data also face accuracy problems. However, there is no good method for the spatial interpolation of precipitation, and research in this area needs to be strengthened.

(4) In the process of zoning in this paper, in addition to the factor selection, there were also problems in calculating the weight of each factor. Because there is no dependent variable, linear regression and other statistical methods cannot be used. The analytic hierarchy process method was used in this article to determine the weights, which are mainly based on the relative importance of each factor and are made by human judgment. Therefore, there are human factors that affect the impact. The results of the zoning will have a certain impact. The analytic hierarchy process is an analytical and discriminative decision-making method, proposed by the American scholar, Saaty, in 1997 [37]. However, although the analytic hierarchy process involves decision-making on data and information, it retains the judgment of the human perception of things and can be completed with less information. The advantage is judging the content of decision-making under the circumstances; however, the disadvantage is that the final result is very dependent on the decisionmaker's judgment of the objective conditions, which has a more obvious subjective intention. When the decisionmaker's perception of things is not accurate enough, it is easy for them to have bias. Therefore, some scholars have proposed methods, such as the multivariate instability index analysis method and the combination weighting method, for weight assignment. We will also conduct further related research.

(5) The current study, based on the meteorological data from 1983 to 2015, reveals the significant reduction in the snow-cover period in Heilongjiang Province over the past 33 years. This reduction is consistent with findings from an analysis of the spatiotemporal variation characteristics of snow in China from 1961 to 2012 [38]. The decrease in the snow-cover period, with interannual variations, has been caused by climate warming over the past 100 years [39]. In addition, the declining trend of the snow period with the latitude variation was similar to that identified in a previous study on the spatiotemporal distribution of snow cover in the northeast region [40]. The snowfall intensity in the east of Heilongjiang Province was higher than in the west, which was consistent with the research results of Huang et al. [30]. The total numbers of blizzards and the total amounts of snowfall in Heilongjiang Province were higher in the east and north, which is similar to the research findings of Zhang [13].