Snow Load Shape Coefficients and Snow Prevention Method for Stepped Flat Roofs

Abstract

Excessive snow load and nonuniform snow deposition are the main factors leading to building collapses. The snow load shape coefficient represents the dimensionless snow load, and its value is related to the unbalanced distribution of snow. The snow load shape coefficients for stepped flat roofs vary greatly in the codes of different regions, which always leads to underestimation of snow loads. We need a widely used standard for snow load shape coefficients. Therefore, through a combination of field measurements and numerical simulations, this study probes the snow accumulation processes and snow load shape coefficients on stepped flat roofs and proposes an equation to calculate snow load shape coefficients and the optimal slope of snow protection for lower roofs. It is found that the maximum snow load shape coefficient emerges at the roof junction with a value of 3.44. The nonuniform length of the snow accumulation is equal to two times the level difference. Based on these, the equation of the snow load shape coefficients is summarized, which is combined with the discrepancies between different codes and the regularity of snow distributions. In this study, the dynamic grid technology under the Eulerian framework is used to successfully predict snow accumulation on stepped flat roofs, and it is noted that snow erosion and deposition are closely related to the location and size of vortexes. Finally, we consider that the ideal slope for the lower roof to prevent snow should be 11°.

1. Introduction

In recent years, extreme snow and ice disasters have occurred frequently around the world, resulting in a significant increase in building damage and collapse [1,2,3]. For example, in 2008, a major snowstorm in southern China caused the collapse of 485,000 buildings in total and a direct economic loss of USD 21 billion. In 2014, a sudden snowfall at the end of winter in China’s Xinjiang province toppled more than 200 houses and damaged 2600 others. Through a large number of post-disaster investigations, it was found that the main reasons for such accidents were excessive snow load and nonuniform snow deposition on the roof [4,5,6]. Unfortunately, there are discrepancies and distinctions between the snow load sections in the codes of different countries. Therefore, it is essential to clarify the distinctions between different national load codes and to conduct intensive research on the snowdrift characteristics on roofs.

The stepped flat roof is a typical snow-sensitive structure. In the current load codes/standards of countries (hereinafter called codes), the snow load values of stepped flat roofs are mainly measured via the field measurements of full-size roofs [7,8,9,10,11]. Since 1956, the Division of Building Research of the National Research Council of Canada has carried out field measurements of snow loads on multiple structures with actual stepped flat roofs, whose results have become the main basis for the modification of snow load codes [11]. By establishing a snowdrift database to statistically explore stepped flat roofs, the USA revised the section on the snow load shape coefficient of stepped flat roofs [10].

The field measurement of the full-size roof is the most intuitive and accurate. However, due to the actual roof’s large span and high snow storage, the final snow pattern is hard to obtain. As a result, reduced-scale field research has gradually developed [12,13,14,15]. Compared with prototype measurement, scaled measurement requires less time and provides more valuable information. Nevertheless, the size effect of snow particles is not yet clear. Fortunately, investigations have shown that equilibrium drifts of models at a scale ratio greater than 1:30 are approximately proportional to the ones of prototype models [13,14,16], and provide suggestions for the snow load shape coefficients of stepped flat roofs [12,16].

Due to the impact of the wind, a lower roof has a more easily unbalanced snow load, resulting in structurally weak areas. Therefore, different approaches are used to reduce the snow load on stepped flat roofs, such as manual operations for snow removal, roof heating systems, and improved roof structures. Manual operations for snow removal are only applicable to the roofs that one can set foot on. In addition, the cost of setting up roof heating systems is high, making them unsuitable for large-scale expansion. In contrast, improving the roof structure is relatively practical. The Guide to the Snow Load Provisions of ASCE 7-16 provides two suggestions for enhancing the structure of stepped flat roofs [17], as shown in Figure 1. The principle of both methods is to add components on the roof to reduce the space available for snowpack. Thus, it is necessary to clarify the snow accumulation processes and the final snow distribution on stepped flat roofs, to provide references and suggestions for the improvement of roof structures.

In conclusion, the snow load codes in different regions are difficult to reconcile, and there is no uniform standard for the stepped flat roof. These codes often obtain different snow load shape coefficients for the same case. Therefore, a harmonized equation to calculate the snow load shape coefficients of stepped flat roofs needs to be proposed. In addition, while raising the slope is well recognized in the scientific field, the optimal slope for the lower roof has still not been specified.

To investigate snow load shape coefficients and the optimal slope for stepped flat roofs, a series of field measurements and numerical simulations were carried out. First, scaled field observations were used to determine the snow pattern on stepped flat roofs. Based on the determination of the snow pattern, the snow load shape coefficients of stepped flat roofs under the final snow status were calculated. The polynomial of the snow load shape coefficients on stepped flat roofs was introduced, which applies to the snow load codes in different countries. Then, dynamic mesh technology modeled with the Eulerian–Eulerian method was successfully used to simulate the snow accumulation processes on stepped flat roofs. Finally, the optimal slope of the lower roof was clarified by comparison to prevent snow accumulation.

2. Field Measurement

2.1. Field Site and Conditions

The field measurement site was in the Mayitas region of Xinjiang, which is located between two mountains 1500 m to 2000 m high. With the prevailing easterly winds, as shown in Figure 2, the maximum instantaneous wind speed could reach 40 m/s due to the “Venturi effect”. In 2018, the annual amount of snowfall in Mayitas was 122.6 mm, which was converted into about 1530 mm of snow depth. All field measurements were conducted at one site to avoid the influence of geographic factors, as shown in Figure 3.

The measurements were carried out on 25 December 2019. The YGY-QXY portable weather station was used to record the daily temperature, humidity, and instantaneous wind speed at a height of 1 m, with a time interval of 1 min. The average relative humidity was 69.2%. The average temperature was −9.6 °C, and the average wind speed was 4.82 m/s. The time history of the wind speed is shown in Figure 4. It is considered that the wind speed conforms to the logarithmic rate:

$$U(z)=\frac{u_{}^{*}}{\kappa }\ln \left(\frac{z}{z_0}\right),$$

(1)

where u* is the friction velocity, κ is the Karman constant equal to 0.42, and z is the vertical height above the ground.

The snow surface roughness length z₀ is calculated using the empirical equation proposed by Pomeroy and Gray (1990) [18], referring to the field measurements:

$$z_0=0.12\frac{u_{*}^2}{2g}.$$

(2)

The friction velocity (Figure 4) is obtained using Equations (1) and (2), and the threshold friction velocity is 0.2 m/s. The atmospheric friction velocity is 0.26 m/s, and the roughness length of the snow surface is 0.4 mm. The measurement lasted for 218 min. Artificial snowfall and natural blowing snow were combined to generate the final pattern of snow distribution on stepped flat roofs.

2.2. Model Parameter

Because of the limitations of the test site and conditions, a model of a stepped flat roof was made of plywood with a scale ratio of 1:25, as shown in Figure 5a. Figure 5b shows the model outline. The dimensions were as follows: 100 mm high (H₂), 300 mm long (L₁), and 600 mm wide at the lower roof and 100 mm high (H₂), 200 mm long (L₂), and 600 mm deep at the upper roof. The length ratio of the lower roof and the upper roof was 3:2, and the height ratio was 1:2.

In order to facilitate the measurement and leveling, a backing board was added to the bottom of the model. In the process of measurement, the status of the snowpack on the roof was maintained. When the snow reached the final distribution, a steel ruler and a drawing square were used to estimate the snow on the center line of the roof. The interval of the measured points was 2 cm. After one measurement, all the snow was removed from the roof and the board.

3. Numerical Simulation

3.1. Simulation Scheme

In this study, the Euler–Euler method was adopted in which both the air phase and snow phase are regarded as a continuum. The proportions of the two phases in the mixture multiphase model are defined by the volume fraction. The continuity equation and momentum equation of the mixed phase are solved by the mixture model.

The continuity equation of the mixed phase is Equation (3):

$$\frac{\partial }{\partial t}\left({\rho }_m\right)+\nabla \cdot \left({\rho }_m{\overrightarrow{v}}_m\right)=0.$$

(3)

The mixed phase density (${\rho }_m$) and the average mass velocity (${\overrightarrow{v}}_m$) are determined using Equations (4) and (5):

$${\rho }_m={\displaystyle \sum }_{k=1}^2{\alpha }_k{\rho }_k,$$

(4)

$${\overrightarrow{v}}_m=\frac{{\displaystyle \sum _{k=1}^2{\alpha }_k{\rho }_k{\overrightarrow{v}}_k}}{{\rho }_m},$$

(5)

where k = 1, 2 represents the principal term (air phase) and secondary term (snow phase), ${\alpha }_k$ is the volume fraction of the kth term, ${\rho }_k$ is the density of the kth term, and ${\overrightarrow{v}}_k$ is the velocity of the kth term.

The momentum equation for the mixture model is expressed as Equation (6):

$$\frac{\partial }{\partial t}\left({\rho }_m{\overrightarrow{v}}_m\right)+\nabla \cdot \left({\rho }_m{\overrightarrow{v}}_m{\overrightarrow{v}}_m\right)=-\nabla p+\nabla \cdot \left[{\mu }_m\left(\nabla {\overrightarrow{v}}_m+\nabla {\overrightarrow{v}}_m^T\right)\right]+{\rho }_m\overrightarrow{g}+\overrightarrow{F}-\nabla \cdot \left({\displaystyle \sum }_{k=1}^2{\alpha }_k{\rho }_k{\overrightarrow{v}}_{dr,k}{\overrightarrow{v}}_{dr,k}\right),$$

(6)

where F is the volume force, and the mixed viscosity and the drift velocity of the kth term are calculated using Equations (7) and (8):

$${\mu }_m={\displaystyle \sum }_{k=1}^2{\alpha }_k{\mu }_k,$$

(7)

$${\overrightarrow{v}}_{dr,k}={\overrightarrow{v}}_k-{\overrightarrow{v}}_m.$$

(8)

The air phase is solved using the Reynolds equation [19]. The turbulence model adopts the realizable k-ε model [20,21], which improves the dissipation rate transport equation compared with the standard k-ε model. Zhao et al. (2016) [22] successfully simulated the wind field and snow distribution on stepped flat roofs with the help of the realizable k-ε model but did not consider the reverse impact of the snowpack on the wind field.

The management of the snow transport and movement in the computational domain is achieved using the snow phase volume fraction (Equation (9)) as well as the erosion and deposition models (Equations (10) and (11)) [23,24,25]. The governing equations are as follows:

$$\frac{\partial }{\partial t}\left({\alpha }_2{\rho }_2\right)+\nabla \cdot \left({\alpha }_2{\rho }_2{\overrightarrow{v}}_m\right)=-\nabla \cdot \left({\alpha }_2{\rho }_2{\overrightarrow{v}}_{\mathrm{dr},2}\right),$$

(9)

$$\begin{array}{cc}{q}_{ero}=A_{ero}\left(u_{*}^2-u_{*t}^2\right)& u_{*}^{}\ge u_{*t},\end{array}$$

(10)

$$\begin{array}{cc}{q}_{{}_{dep} }=\varphi w_f\frac{u_{*t}^2-u_{*}^2}{u_{*t}^2}& u_{*}^{}<u_{*t},\end{array}$$

(11)

where A_ero is the erosion constant, equal to 7.0 × 10⁻⁴ [23]; u_* is the wall friction velocity; u_*t is the threshold friction velocity of snow drifting; and Φ is the snow mass concentration.

3.2. Computational Parameters

Table 1. Parameters of the numerical simulation and field measurements.

Parameter	Field Measurement	Numerical Simulation
Air density ρ (kg/m3)	1.225	1.225
Snow density ρs (kg/m3)	120~300	150
Snow particle diameter ds (μm)	165 (Expected value)	165
Falling speed wf (m/s)	0.2~0.5	0.2 (h > 0.1 m), 0.35 (h ≤ 0.1 m)
Threshold friction velocity u*t (m/s)	0.2	0.2

presents the major parameters in the field measurements and the numerical simulation. Ma et al. (2021b) [26] found that the probability distribution of the equivalent diameter of blowing snow particles meets the gamma distribution Ga(α, β), where α and β are the shape parameter and scale parameter (α = 2.94, β = 100.35), respectively, and the expected integrated snow particle size is 165 μm [27]. The measured snow density ranges from 120 kg/m³ to 300 kg/m³, and the snow density was 150 kg/m³ in the numerical simulation [28,29,30]. During snow transport, the diameter and the falling velocity of snow particles tend to be smaller in the higher part above the snow surface. Therefore, a lower falling speed (w_f = 0.2 m/s) [28,30,31,32] was adopted in the higher part (h > 0.1 m), while in the lower part (h ≤ 0.1 m), a higher falling speed (w_f = 0.35 m/s) was used [33,34]. The threshold friction velocity in the numerical simulation was 0.2 m/s [28,30], as illustrated in Figure 4.

3.3. Calculation Settings

The study was carried out using ANSYS. In this work, a 2D model was adopted for simulation. The size of the numerical simulation model was the same as that of the model used in the field measurements, and the dimensions of the computational domain were 8000 mm (16 L) × 3000 mm (15 H₂), as shown in Figure 6. The length of the rear flow field was 5000 mm (10 L) to fully develop the snow flow. In the processes of snow drift, the snow concentration above 1 m from the snow surface was almost 0; so, the height of the calculation domain was considered to be 3000 mm (15 H₂).

The number of core grids accounts for more than 50% of the total grid. In order to maximize the quality of the grid around the model, an even grid with a height of 5 mm was adopted around the model. The mesh growth factor was 1.1 in the region far from the model, and the y⁺ value was controlled at about 80. The total number of meshes was about 30,000, the quality and orthogonality of the meshes were about 1, and the maximum aspect ratio of the meshes was less than 13.

The measured logarithmic wind speed profile was imposed by Equation (1) at the inlet, and the turbulent kinetic energy and dissipation rate were obtained using Equations (12) and (13) [35], with C_μ, the empirical constant, equal to 0.09.

$$k\left(z\right)=\frac{u_{ABL}^{*2}}{\sqrt{C_{\mu }}}$$

(12)

$$\epsilon \left(z\right)=\frac{u_{ABL}^{*3}}{\kappa z}$$

(13)

The inlet snow phase was defined by the volume fraction, which refers to Equations (14) and (15) derived by Kang et al. (2018) [34]:

$$f_{sal }=\frac{\rho }{3.29{\rho }_s{u}_{*}}\left(1-\frac{u_{*t}^2}{u_{*}^2}\right),$$

(14)

$$f_{susp }=f_{sal }{\left(\frac{z}{h_{sal }}\right)}^{-\left[w_f/\left(\kappa u_{*}\right)\right]},$$

(15)

where h_sal is the height of the saltation layer, which can be calculated using Equation (16) [18]:

$$h_{sal}=1.6\frac{u_{*}^2}{2g}.$$

(16)

The boundary of the ground and the model are set as the no-slip wall boundary, and the snow surface contour is updated in real time by using dynamic grid technology based on the transient method. The dynamic grid technology can consider the effect of the snow boundary on the airflow near the model and the trajectory of subsequent snow particles, to further improve the simulation accuracy [33,34]. The principle of the dynamic grid is shown in Figure 7.

For the values of the wall roughness height (K_S) and roughness constant (C_S), this work set the roughness constant C_S = 2.0 [34]. Based on Equation (17) proposed by Blocken et al. (2007) [36], the roughness height of the wall was obtained as K_S = 2 mm, which satisfies K_S < y_p, and y_p is the center height of the minimum grid near the ground (y_p = 2.5 mm).

$$K_S=\frac{9.793z_0}{C_S}$$

(17)

The upper boundary of the domain adopts a symmetry boundary, and the outlet is set as the pressure outlet. The calculation settings are summarized in

Table 2. Calculation settings of the numerical simulation.

Calculation Settings
Computational domain	8000 mm (16 L) × 3000 mm (15 H2)
Mesh division	Minimum grid height: 5 mm. Grid growth rate: 1.1. Total grid: 30,000. y+: 80. Mesh quality: the quality and orthogonality are about 1, and the aspect ratio is less than 13.
Inlet boundary	Air phase: Equation (1) and Equations (12) and (13). Snow phase: Equations (14)~(16).
Ground and model boundary	Wall: CS = 2.0, KS = 2 mm.
Outlet boundary	Pressure outlet
Upper boundary	Symmetry

.

To increase the credibility of our work, the methodological flow diagram is shown in Figure 8.

4. Results and Discussion

4.1. The Final Pattern of Snow Accumulation

The corresponding wind speed and temperature in the field measurements are illustrated in Figure 9. When the lower roof was on the windward side, the mean wind speed was 4.68 m/s, and the mean temperature was −10.3 °C. It took about 13 min to reach the final pattern of snow distribution. On the other hand, the windward snowdrift reached its final pattern in about 18 min, during which the mean wind speed was 4.91 m/s, and the mean temperature was −9.3 °C. The field measurements were conducted from 13:07 to 14:13, and the wind speed and temperature were relatively stable.

Figure 10 presents the measured snow deposition on the roof. In Figure 11, the horizontal axis denotes the distance from the upper wall x normalized by the level difference H, and the vertical axis denotes the dimensionless height h/H. The origin of the horizontal axis is located at the junction of the upper roof and the lower roof.

As indicated in Figure 10a and Figure 11a, when the lower roof was on the windward side, there was almost no snow accumulation on the upper roof. The final snow deposition on the lower roof was similar to a triangle, with the peak snow depth located at the roof junction (x/H = 0) and aligned with the outer edge of the upper roof (h/H = 2.0). The snow accumulation leveled off at x/H = 1.4. Due to the wind separation and acceleration occurring at x/H = 3, the lower roof presented a snow-free condition at x/H = 2.95~3.

Figure 10b and Figure 11b show the snow distribution when the upper roof was on the windward side. There was a small amount of snow on the upper roof. The peak snow depth was about 1/2H at x/H = 0. Similar to the appearance in Figure 10a and Figure 11a, the wind accelerated at the windward corner of the lower roof (x/H = −2), and the snow-free area occurred at x/H = −2~−1.8. The snow distribution on the lower roof was approximately uniform, and the dimensionless snow depth was around 0.4 at x/H = 0.2~2.4.

4.2. Snow Load Shape Coefficients

The snow depth d on the roof is normalized by the variation in snow depth Δd on the ground during the measurement. The dimensionless values obtained are regarded as snow load shape coefficients (μ), as shown in Figure 12.

The windward snow load shape coefficients on the lower roof gradually decreased with the increase in x/H. The snow load on the roof at x/H = 1.6 was the same as the ground snow load. The maximum snow load shape coefficient occurred at x/H = 0, μ_max = 3.44. As indicated in Figure 12a, the snow load shape coefficients demonstrated a strong regularity, which was fitted with a quadratic polynomial of x/H to obtain Equation (18) (R² = 0.9991):

$$\mu =0.3346{\left(x/\mathrm{H}\right)}^2-2.0547\left(x/\mathrm{H}\right)+3.4426.$$

(18)

The snow load shape coefficient of the lower roof when the upper roof was on the windward side (Figure 12b) was about 1 at x/H = 0.2~2.4. The snow load shape coefficient gradually decreased to 0.25 at x/H ≥ 2.4. However, this area was relatively small compared with the overall low roof; so, the snow load shape coefficient was still unified as 1. The maximum value of the shape coefficient on the lower roof occurred at the junction of the roofs (x/H = 0), μ_max = 1.25. The snow load shape coefficient on the upper roof was smaller than that on the lower roof, and the maximum value was about 0.5.

The snow load on stepped flat roofs is specified in detail in many countries. However, the definitions of the snow load shape coefficient are different. Comparing the codes of different countries, the formula of the snow load values of roofs can be expressed in Equation (19) (the ASCE can also be derived to obtain this form), in which S is the snow load on the roof, μ is the snow load shape coefficient, and S₀ is the ground snow load.

$$S=\mu S_0$$

(19)

In this paper, the different expressions of snow load in different national load codes are unified into Equation (19). As listed in

Table 3. Formulas of snow load shape coefficients in the codes of different countries.

Code	Calculation Formula of μ	Parameter Description
GB50009 2012 [9]	$\begin{array}{c}\begin{array}{cc}\mu =\left(L_1+L_2\right)/2H& (2.0\le \mu \le 4.0)\end{array}\\ \begin{array}{cc}{x}_d=2H& (4\mathrm{m}\le x_d\le 8\mathrm{m})\end{array}\end{array}$	L1 is the lower roof span. L2 is the upper roof span. H is the height difference between the stepped flat roofs. xd is the length of the nonuniform snow distribution. μb is the snow uniform distribution coefficient. μw is the drift coefficient. μs is the slide coefficient. Ca and μi are the shape coefficients of the snow load. Ce is the exposure coefficient. Ct is the temperature coefficient. Cs is the roof slope coefficient. γ is the weight density of snow. Is is the building importance factor. lc is the characteristic length of stepped flat roofs. hb is the uniform snow thickness. hd is the snowdrift thickness. hc is the height difference from the uniform snow surface to the upper roof. Sb is the uniformed snow load.
ASCE 2017 [10]	$\begin{array}{c}\begin{array}{c}\mu =\frac{\gamma \left(h_d+h_b\right)}{S_0}\\ \begin{array}{cc}\frac{h_d}{\sqrt{I_s}}=\left(0.43\sqrt[3]{L_1}\sqrt[4]{S_0+10}\right)-1.5& \left(h_d\le h_c\right)\end{array}\end{array}\\ h_b=S_b/\gamma \\ \begin{array}{cc}\gamma =0.13S_0+14& \left(\gamma \le 30\mathrm{lb}/{\mathrm{ft}}^3\right)\end{array}\\ \begin{array}{cc}{x}_d=\left\{\begin{array}{c}\begin{array}{cc}4h_d& \left(h_d\le h_c\right)\end{array}\\ \begin{array}{cc}4h_d^2/h_c& \left(h_d>h_c\right)\end{array}\end{array}\right.& \left(x_d\le 8h_c\right)\end{array}\end{array}$
EN 1991-1-3 2003 [8]	$\begin{array}{c}\begin{array}{l}\mu ={\mu }_i{C}_{\mathrm{e}}{C}_t\\ {\mu }_2={\mu }_{\mathrm{w}}+{\mu }_{\mathrm{s}}\end{array}\\ \begin{array}{cc}{\mu }_{\mathrm{w}}=\left(b_1+b_2\right)/2h& \left({\mu }_{\mathrm{w}}\le \gamma H/S_0\right)\end{array}\\ {\mu }_{\mathrm{s}}=\left\{\begin{array}{ll}0\hfill & \alpha <{15}^{\circ }\hfill \\ 0.5{\mu }_1\hfill & \alpha \ge {15}^{\circ }\hfill \end{array}\right.\\ \begin{array}{cc}{x}_d=2H& \left(5\mathrm{m}\le x_d\le 15\mathrm{m}\right)\end{array}\end{array}$
NBCC 2015 [11]	$\begin{array}{c}\mu =I_s{\mu }_b{C}_e{C}_s{C}_a\\ {\mu }_b=\left\{\begin{array}{c}\begin{array}{cc}0.8& l_c\le \left(\frac{70}{C_e^2}\right)\end{array}\\ \begin{array}{cc}\frac{1}{C_e}\left[1-\left(1-0.8C_e\right)\exp \left(-\frac{l_c{C}_e^2-70}{100}\right)\right]& l_c>\left(\frac{70}{C_e^2}\right)\end{array}\end{array}\right.\\ C_a=\left\{\begin{array}{c}\begin{array}{cc}{C}_{a0}-\left(C_{a0}-1\right)\left(x/x_d\right)& 0\le x\le x_d\end{array}\\ \begin{array}{cc}1.0& x>x_d\end{array}\end{array}\right.\\ \begin{array}{cc}{x}_d=5\frac{{\mu }_b{S}_0}{\gamma }\left(C_{a0}-1\right)& C_{a0}=\min \left(\beta \frac{\gamma H}{{\mu }_b{S}_0},\frac{F}{{\mu }_b}\right)\end{array}\end{array}$
AIJ 2004 [7]	$\begin{array}{c}\mu ={\mu }_b+{\mu }_w+{\mu }_s\\ {\mu }_s=\left\{\begin{array}{ll}0\hfill & \alpha <{10}^{\circ }\hfill \\ {\mu }_b\hfill & \alpha \ge {25}^{\circ }\hfill \end{array}\right.\end{array}$

, the set of dimensionless coefficients that are multiplied by the S₀ in each code is defined as the snow load shape coefficient. For easier comparison, the parameters with the same meanings in the formulas of different countries have been standardized.

The greater values of snow load shape coefficients at the same dimensionless distance are selected in two measurements, which are used to obtain the expressions for the snow load shape coefficients in this work as Equation (20):

$$\mu =\left\{\begin{array}{cc}0.3346{\left(x/\mathrm{H}\right)}^2-2.0547\left(x/\mathrm{H}\right)+3.4426& 0\le x/\mathrm{H}\le 1.6\\ 1& x/\mathrm{H}>1.6\end{array}\right..$$

(20)

The fitting equation is compared with the snow load shape coefficients in different codes in Figure 13. Each code has fewer definitions of snow distribution on the upper roof, and the measured snow depth of the lower roof is much greater than that of the upper roof. Moreover, the snow load shape coefficients on the upper roof are far less than 1. Therefore, only the snow load shape coefficients for the lower roof are compared.

For stepped flat roofs, different national load codes hold that the maximum snow load shape coefficient emerges at the roof junction (x = 0), which is consistent with the measured results in this paper. However, the codes of different countries assume that the nonuniform distribution of snow load is a linear regularity. The measurement results justify that the law of quadratic polynomials is more suitable for expressing the snow distribution coefficients. Additionally, the maximum snow load shape coefficient in this paper is higher than that of other codes. The snow load shape coefficients of the AIJ code are the lowest in the snow load codes, with a maximum of only about 1. The ASCE code states that the unbalanced length of the snowpack only accounts for 1/6 of the roof span with snow load shape coefficients of less than 2. The snow load shape coefficients in NBCC are the highest among the other codes, whose maximum value is 2.7, but the nonuniform length of the snow accumulation is equal to the level difference (H). The EN and GB codes all consider that the length for the nonuniform area is 2/3 (2H) of the roof span, and the maximum snow load shape coefficient in the two codes is 2.5. The maximum snow load shape coefficients and the unbalanced length of snow are summarized in

Table 4. Comparison of the maximum snow load shape coefficients and the unbalanced length of snow.

Code or Equation	Maximum μ	Dimensionless Unbalanced Length (x/H)
AIJ	1	2
ASCE	2	0.5
NBCC	2.7	1
EN	2.5	2
GB	2.5	2
Equation (20)	3.44	1.6
Equation (21)	3.44	2

.

It is evident that the values of the snow load shape coefficients in different codes vary significantly due to the differences in the environmental and geographical characteristics of each country. This poses challenges in identifying the snow load of stepped flat roofs in practical projects. Therefore, the maximum snow load shape coefficients at each position and the maximum nonuniform length of deposition are summarized in Equation (21), which refers to the codes of different countries, and Equation (20).

$$\mu =\left\{\begin{array}{cc}-1.22\left(x/\mathrm{H}\right)+3.44& 0\le x/\mathrm{H}\le 2\\ 1& x/\mathrm{H}>2\end{array}\right.$$

(21)

Equation (20) states that the unbalanced length of the snow load on the lower roof is twice the level difference. The values of the snow load shape coefficient conform to linear regularity, with the maximum located at x = 0, and the maximum is about 3.44.

4.3. Snow Accumulation Processes for Stepped Flat Roofs

Based on the correspondences between the measurement results and codes, differences can be noted in the snow load shape coefficient values on stepped flat roofs. To explore the reasons, Figure 14 compares the snow accumulation processes on the stepped flat roof and flow fields around the roof at different times. The simulation time is 780 s, which is the same as the field-measured time. The dimensionless length of the erosion zone is defined as L_e, and the time t of the snow accumulation is nondimensionalized to T = t/780. Correspondingly, T = 1 represents the time when the snow evolved into its final form. The snow shape on the lower roof was almost rectangular at T₁ = 0.064, and the erosion and the deposition of snow were not obvious at this time. However, via the flow field analysis, it can be observed that a vortex with a dimensionless length of about 2 was formed on the lower roof, from which it can be inferred that this area was the erosion region. When the snow developed to T₂ = 0.128, obvious erosion was observed due to the vortex on the lower roof. The dimensionless range of the erosion length was about 0~1. Over time, the erosion region gradually shrank. The vortex scale and the erosion length decreased to 0.5 at T₃ = 0.256. The reverse flow on the lower roof disappeared at T₄ = 0.385 when erosion was also not observed. It was concluded that erosion could be observed at T₂ = 2T₁, and the erosion length was reduced by half at T₃ = 4T₁. This indicates a linear relationship between the erosion length L_e and the time T of the snow accumulation on stepped flat roofs; that is, L_e = CT, where C is the proportional coefficient. The snow accumulation reached the final pattern at T₅ = 1. In the processes of snow accumulation, a flow separation preventing deposition occurred at the end of the lower roof, which verifies the results observed in the measurement.

The simulation results of snow development are compared with the measurement results [12,14,22], as shown in Figure 15. The dotted lines indicate the numerical simulation results at different dimensionless times, while the solid lines represent the measurement results. The relative results are marked with the same color.

Figure 15 indicates that the simulation results are consistent with the field measurement results, and the erosion area and the peak depth of snow are reflected in the simulation results. Further comparison illustrates that different snow distributions in the literature correspond to different moments of snow accumulation on stepped flat roofs. Moreover, the corresponding times are all less than one-third of the snow accumulation time. This also proves that it is difficult for the snow on prototype structures to reach a stable state due to the influence of the measurement time and the environment.

4.4. Snow Prevention Method for Stepped Flat Roofs

According to the results of the measurements and numerical simulation, the final shape of the windward snow accumulation on the lower roof is approximated as a triangle. It can be predicted that increasing the angle of the lower roof can reduce the snow depth. Therefore, when the angle of the lower roof β is 1°~14°, the final snow distribution of the lower roof is obtained as shown in Figure 16. As the angle of the lower roof increases, the snow profile on the lower roof is virtually unchanged, but the amount of snow accumulation significantly decreases. To quantitatively analyze the change in snow accumulation, the snow removal rate σ is defined as Equation (22):

$$\sigma =\frac{A_0-A_{\beta }}{A_0},$$

(22)

where A_β is the snow accumulation on the lower roof at different angles, and A₀ represents the initial amount of snow on the lower roof.

The snow removal rate σ and the maximum snow depth D on the lower roof are depicted in Figure 17. When β = 0°~3°, the amount of snow accumulation on the lower roof decreases by more than 10%, and the snow depth significantly decreases with the reduction in the snow accumulation. Subsequently, the snow removal rate and snow depth show a linear trend in the range of β = 3°~11°. When β > 11°, as the snow removal rate increases, the snow depth decreases slowly. The snow removal on the lower roof is more than 90% at β = 12°, and the maximum dimensionless snow depth is less than 0.2. When β increases from 13° to 14°, the snow accumulation on the lower roof reduces by only 1.8%, with the snow depth less than 0.1H.

As the angle gradually increases, the amount of snow accumulation on the roof decreases. However, the increase in the slope angle will extend the roof, which results in higher construction costs and resource waste. Therefore, the snow load shape coefficients of the lower roof are calculated as shown in Figure 18. The maximum snow load shape coefficient of the lower roof gradually decreases with the increase in the roof slope angle. When the slope increases by 1°, the maximum snow load shape coefficient decreases by about 0.2. The snow load on the lower roof is less than the ground snow load at β = 11°. When β ≥ 13°, the maximum snow load shape coefficient of the lower roof has reached below 0.5.

Finally, compositing the snow removal rate σ, maximum snow depth D, and snow load shape coefficient μ on the lower roof, the slope of the lower roof is suggested to be 11°. When β = 11°, the snow removal rate on the lower roof is more than 85%, the maximum snow depth is about 0.25H, and the maximum snow load shape coefficient is about 0.85.

5. Conclusions

This study combines field measurements and numerical simulations to study the snow accumulation processes and final patterns of snow distribution on a stepped flat roof. By integrating the measured results with different codes, a formula for the snow load shape coefficient is presented. Moreover, a viable design option to minimize the potential for any drift is proposed for snow prevention on stepped flat roofs. Through these efforts, the following conclusions are obtained:

(1)

During the snow accumulation processes on stepped flat roofs, an erosion region whose length varies linearly with time is formed on the lower roof due to reverse flow. The erosion disappears when the snow accumulation process is almost half-finished. As the snow develops to its final pattern, no measurable snow is observed on the upper roof, while the snow shapes on the lower roof vary with different wind directions. The final windward snow on the lower roof is distributed in a triangle. When the upper roof is on the windward side, the final snow distribution is close to uniform.

(2)

The codes of different countries have various snow load shape coefficients on stepped flat roofs, but it is generally accepted that the nonuniform snow distribution on stepped flat roofs is close to linear. The measured snow load shape coefficients are more in line with a quadratic function. Through the combination of measured results and different codes, it is recommended that the value of the snow load shape coefficient conform to a linear regularity: $\mu =\left\{\begin{array}{cc}-1.22\left(x/\mathrm{H}\right)+3.44& 0\le x/\mathrm{H}\le 2\\ 1& x/\mathrm{H}>2\end{array},\right.$ where x/H is the distance normalized by the level difference H.

(3)

An increase in the slope can effectively reduce the snow accumulation on stepped flat roofs. As the slope increases, the snow removal rate of the roof gradually increases, while the maximum snow depth and snow load shape coefficients subsequently decrease. When the snow load on the roof is less than the ground snow load, there is an optimal slope of the lower roof, making the snow removal rate the maximum and the area of the lower roof the minimum. Finally, considering snow removal efficiency and economic factors, the slope of the lower roof is recommended to be 11°.