METHOD FOR WHOLE-PROCESS NUMERICAL SIMULATION AND HAZARD FORECAST OF MOUNTAIN DISASTER

Abstract

A method for a whole-process numerical simulation and hazard forecast of a mountain disaster is provided. The method includes: S1, a high space-time rainfall forecast of a mountain area; S2, a hydrodynamic process and numerical simulation: establishing a hydrodynamic process model and solving the hydrodynamic process model; S3, a motion model and numerical simulation of a mountain torrent and debris flow disaster; and S4, a risk analysis and hazard forecast of a small watershed disaster. The present invention predicts disaster hazard and dynamically and quantitatively evaluates risk loss according to a whole-process scenario simulation of the disaster driven by a climate forecast result, improves current disaster level forecasts to hazard forecasts, and serves for accurate disaster preventions and accurate rescues.

Claims

1. A method for a whole-process numerical simulation and a hazard forecast of a mountain disaster, comprising: S1, a high space-time rainfall forecast of a mountain area: interpolating forecast data having a space resolution of 9 KM to 0.01°*0.01° by means of a bilinear interpolation method, and establishing an empirical relation between a physical quantity and an elevation according to terrains to describe an influence of the terrains on a precipitation; S2, a hydrodynamic process and numerical simulation: establishing a hydrodynamic process model and solving the hydrodynamic process model; S3: a motion model and numerical simulation of a mountain torrent and debris flow disaster, wherein a depth-averaged continuous medium equation for a facies averaging rainfall-induced debris flow is expressed as: $\frac{\partial h}{\partial t}+\frac{\partial \mathrm{hu}}{\partial x}+\frac{\partial \mathrm{hv}}{\partial y}=R-I+\frac{E}{1-p}$ $\frac{\partial \mathrm{hu}}{\partial x}+\frac{\partial \left({\mathrm{hu}}^2+0.5{\mathrm{gh}}^2\right)}{\partial x}+\frac{\partial \mathrm{huv}}{\partial y}=-\mathrm{gh}\frac{\partial z_b}{\partial x}-\frac{{\tau }_{\mathrm{fx}}}{\rho }-\frac{\left({\rho }_s-{\rho }_f\right){\mathrm{gh}}^2}{2\rho }\frac{\partial c}{\partial x}+\frac{\left({\rho }_b-\rho \right)\mathrm{uE}}{\rho \left(1-p\right)}$ $\frac{\partial \mathrm{hv}}{\partial t}+\frac{\partial \mathrm{huv}}{\partial x}+\frac{\partial \left({\mathrm{hv}}^2+0.5{\mathrm{gh}}^2\right)}{\partial y}=-\mathrm{gh}\frac{\partial z_b}{\partial y}-\frac{{\tau }_{\mathrm{fy}}}{\rho }-\frac{\left({\rho }_x-{\rho }_f\right){\mathrm{gh}}^2}{2\rho }\frac{\partial c}{\partial y}+\frac{\left({\rho }_b-\rho \right)\mathrm{vE}}{\rho \left(1-p\right)}$ $\frac{\partial \mathrm{hc}}{\partial t}+\frac{\partial \mathrm{huc}}{\partial x}+\frac{\partial \mathrm{hvc}}{\partial y}=E$ $\frac{\partial z_b}{\partial t}=-\frac{E}{1-p}$ wherein t represents time, x and y are horizontal coordinates, h represents a fluid depth, u and v are components of a fluid depth-averaged speed in an x direction and a y direction, respectively, c is a depth-averaged solid phase concentration, g is gravitational acceleration, R is a reduced rainfall intensity after vegetation interception, I is an infiltration rate, Z.sub.b is a substrate elevation, ρ is a density of a solid-liquid mixture, ρ=cρ.sub.s+(1−c)ρ.sub.f, ρ.sub.s and ρ.sub.f are concentrations of a solid phase and a liquid phase respectively, ρ.sub.b is a density of a saturated substrate, ρ.sub.b=(1−d)ρ.sub.s+pρ.sub.f, p is the porosity of a substrate material, τ.sub.fx and τ.sub.fy are substrate resistances in the x direction and the y direction respectively, and E is an erosion rate of the substrate; and S4, a risk analysis and hazard forecast of a small watershed disaster, further comprising: computing a comprehensive disaster-causing risk degree of a disaster: determining the comprehensive disaster-causing risk degree of the disaster according to composite hazard characteristics of impact and burying of the mountain disaster; analyzing a vulnerability: computing a vulnerability of disaster-bearing bodies according to values of the disaster-bearing bodies and vulnerability indexes of the disaster-bearing bodies; and computing a risk degree of the disaster: determining the risk of the disaster on the basis of a numerical simulation result, wherein the risk of the disaster is a comprehensive function of a comprehensive hazard of the disaster, the vulnerability of the disaster-bearing bodies and exposure of the disaster-bearing bodies; and a governing equation for the hydrodynamic process model is: $\frac{\partial h}{\partial t}+\frac{\partial \mathrm{hu}}{\partial x}+\frac{\partial \mathrm{hv}}{\partial y}=R-V-I$ $\frac{\partial \left(1/2{\mathrm{gh}}^2\right)}{\partial x}=-\mathrm{gh}\frac{\partial z_b}{\partial x}+S_{\mathrm{fx}}$ $\frac{\partial \left(1/2{\mathrm{gh}}^2\right)}{\partial y}=-\mathrm{gh}\frac{\partial z_b}{\partial y}+S_{\mathrm{fy}}$ wherein t is time, x and y are the horizontal coordinates, h represents the fluid depth, u and v represent the speed components of the fluid depth-averaged speed in the x direction and the y direction, respectively, R is the reduced rainfall intensity after the vegetation interception, I is the infiltration rate, V is the vegetation interception, g is the gravitational acceleration, Z.sub.b is the substrate elevation, S.sub.fx and S.sub.fy represent substrate friction terms in the x direction and the y direction, respectively, and S.sub.fx and S.sub.fy are expressed by Manning models as: $S_{\mathrm{fx}}=g\frac{n^{2}u\sqrt{u^2+v^2}}{h^{1/3}}$ $S_{\mathrm{fy}}=g\frac{n^{2}v\sqrt{u^2+v^2}}{h^{1/3}}$ wherein n is a Manning coefficient, h represents the fluid depth, u and v represent the speed components of the fluid depth-averaged speed in the x direction and the y direction, respectively, and g is the gravitational acceleration; an Aston vegetation rainfall interception model used by the vegetation interception V is expressed as: $V=S_{\max }\left(1-e^{-k\frac{P_c}{S_{\max }}}\right)$ wherein S.sub.max is maximum interception of vegetation, and according to computation of different types of vegetation, P.sub.c is accumulated rainfall, and k is a parameter related to a canopy density of the vegetation, and is expressed as:
k=1−e.sup.−(Co*LAI), wherein Co is the canopy density of the vegetation, and LAI is a leaf area index; and the infiltration rate I is expressed by a slope saturation infiltration model as: $I=\frac{\mathrm{df}}{\mathrm{dt}}=k_s\left[\left({\psi }_f+h\right)\frac{{\theta }_s-{\theta }_i}{f}+1\right]$ wherein k.sub.s is a saturated hydraulic coefficient, ψ.sub.f is a matrix suction water head at a front end of a wetting front, θ.sub.s is saturated moisture content of soil, θ.sub.i is initial moisture content of the soil, and f is an accumulated infiltration depth.

2. The method according to claim 1, wherein in step S1, the bilinear interpolation method comprises: computing an attribute value at a point P.sub.1=(x, y), carrying out linear interpolation in the x direction to obtain $f\left(R_1\right)\approx \frac{x_2-x}{x_2-x_1}f\left(Q_{11}\right)+\frac{x-x_1}{x_2-x_1}f\left(Q_{21}\right),\mathrm{wherein}$ $R_1=\left(x,1_{}\right),\mathrm{and}$ $f\left(R_2\right)\approx \frac{x_2-x}{x_2-x_1}f\left(Q_{12}\right)+\frac{x-x_1}{x_2-x_1}f\left(Q_{22}\right),\mathrm{wherein}$ $R_2=\left(x,2_{}\right),$  and then carrying out linear interpolation in the y direction to obtain $f\left(P_1\right)\approx \frac{2_{}-}{2_{}-1_{}}f\left(R_1\right)+\frac{-1_{}}{2_{}-1_{}}f\left(R_2\right)$ and to further obtain a desired result f(x, y): $f\left(x,\right)\approx \frac{f\left(Q_{11}\right)}{\left(x_2-x_1\right)\left(2_{}-1_{}\right)}\left(x_2-x\right)\left(2_{}-\right)+\frac{f\left(Q_{21}\right)}{\left(x_2-x_1\right)\left(2_{}-1_{}\right)}\left(x-x_1\right)\left(2_{}-\right)+\frac{f\left(Q_{12}\right)}{\left(x_2-x_1\right)\left(2_{}-1_{}\right)}\left(x_2-x\right)\left(-1_{}\right)+\frac{f\left(Q_{22}\right)}{\left(x_2-x_1\right)\left(2_{}-1_{}\right)}\left(x-x_1\right)\left(-1_{}\right)$ wherein f represents a corresponding attribute value at a certain point, and Q.sub.11, Q.sub.12, Q.sub.21 and Q.sub.22 represent point positions at (x.sub.1, y.sub.1), (x.sub.1, y.sub.2), (x.sub.2, y.sub.1) and (x.sub.2, y.sub.2).

3. The method according to claim 1, wherein in step S2, the establishing of the hydrodynamic process model comprises: carrying out a deep integral simplification on the basis of a Navier-Stokes equation and ignoring convection terms in a momentum equation to obtain a diffusion wave model, quantitatively computing a hydrodynamic process of a small watershed, introducing the Aston vegetation rainfall interception model and a Green-Ampt slope saturation infiltration model on the basis of the diffusion wave model to establish the hydrodynamic process model, and simulating a whole physical event from rainfall to vegetation interception and slope infiltration and then to runoff generation and motion.

4. The method according to claim 1, wherein in step S2, the solving of the hydrodynamic process model comprises: carrying out a numerical solution by means of a first-order windward difference scheme to carry out fine-grained parallelization on a program.

5. The method according to claim 1, wherein in step S4, the comprehensive disaster-causing risk degree of the disaster is computed by the following model: $H=H_e+H_d,$ $H_d=\frac{N_{i,j}\Delta V}{A}$ $H_e=A.Math.\underset{t>0}{\max }\left[\left(u^2+v^2\right)h\rho \right]$ wherein H is the comprehensive disaster-causing risk degree of the disaster, and H.sub.e is a hazard caused by impact damage and is expressed by maximum kinetic energy of a disaster-causing body; H.sub.d is a hazard caused by burying and is expressed by a maximum burying depth of the disaster-causing body; N.sub.i, j is equal to the number of particles in a control grid centered on a point (i, j); ΔV is a volume of the particle; A is a grid area; and u and v represent the speed components of the fluid depth-averaged speed in the x direction and the y direction, respectively, h is the fluid depth, and p is the density of the solid-liquid mixture.

6. The method according to claim 1, wherein in step S4, the computing of the vulnerability of the disaster-bearing bodies according to the values of the disaster-bearing bodies and the vulnerability indexes of the disaster-bearing bodies comprises: describing vulnerability of the i.sup.th type of disaster-bearing bodies as:
V.sub.i=V(u).sub.i×C.sub.i, wherein V.sub.i is the vulnerability of the i.sup.th type of disaster-bearing bodies, V(u).sub.i is a comprehensive value of the i.sup.th type of disaster-bearing bodies, and C.sub.i is a vulnerability index of the i.sup.th type of disaster-bearing bodies, wherein quantification of the comprehensive value V(u) of the disaster-bearing bodies depends on an average unit price D.sub.e of the disaster-bearing bodies and an actual area A.sub.e affected by the disaster, which is expressed as:
V(u)=D.sub.e*A.sub.e, and the vulnerability index C of the disaster-bearing bodies is expressed as: $C=\frac{h}{H_c}$ wherein h is the fluid depth, H.sub.c is a geometric height of the disaster-bearing bodies, and C is valued as 1 when $\frac{h}{H_c}\ge 1.$

7. The method according to claim 1, wherein in step S4, the risk degree of the disaster is computed according to the following equation: $\mathrm{Ra}=f\left(H,V,E\right)=\underset{i=1}{\overset{n}{.Math.}}{H}_i\times V_i\times E_i$ wherein Ra is the risk degree of the disaster and is expressed by a certain numerical value ranging from 0 to 1; H.sub.i is a comprehensive disaster-causing risk degree of the i.sup.th type of disasters and is expressed by a certain numerical value ranging from 0 to 1; V.sub.i is a vulnerability degree of the i.sup.th type of disaster-bearing bodies and is determined by damage degrees and values or the number of the disaster-bearing bodies; and E.sub.i is an exposure degree of the i.sup.th type of disaster-bearing bodies, and is expressed by one of numerical values ranging from 0 to 1.

Description

BRIEF DESCRIPTION OF DRAWINGS

[0051] In order to more clearly illustrate the specific implementations of the present invention or the technical solutions in the prior art, a brief introduction to the accompanying drawings required for the description of the specific implementations or the prior art will be provided below. Like elements or parts are generally identified by like reference numerals in all the drawings. In the drawings, elements or parts are not necessarily drawn to actual scale.

[0052] FIG. 1 is a flow chart of a high-resolution data processing technology in the present invention;

[0053] FIG. 2 is a schematic diagram of bilinear interpolation in the present invention;

[0054] FIG. 3 is a technical route diagram of a hydrodynamic module in the present invention;

[0055] FIG. 4 is a topographic map of a research area in the present invention;

[0056] FIG. 5 is a land utilization type map in the present invention;

[0057] FIG. 6 is a flow process graph in the present invention;

[0058] FIG. 7 is a map of a water depth simulation result of a Longxi river in the present invention;

[0059] FIG. 8 is a vulnerability grading map of a Longxi river watershed; and

[0060] FIG. 9 is a map of a risk grading result of a mountain disaster in the Longxi river watershed.

DETAILED DESCRIPTION OF EMBODIMENTS

[0061] The technical solutions in the embodiments of the present invention are clearly and completely described below. Apparently, the embodiments described are merely some rather than all of the embodiments of the present invention. On the basis of the embodiments in the present invention, all other embodiments acquired by those of ordinary skill in the art without making creative efforts fall within the scope of protection of the present invention.

[0062] It should be noted that in all directional indications (such as top, bottom, left, right, front, rear . . . ) in the embodiments of the present invention are only used for explaining a relative position relation, motion, etc. between components under a certain specific attitude (as shown in the figures), and if the specific attitude is changed, the directional indications will be changed accordingly.

[0063] In addition, descriptions of “first,” “second”, etc. involved in the present invention are merely used for descriptive purposes and are not to be construed as indicating or implying their relative importance or implicitly specifying the number of indicated technical features. Thus, a feature defined with “first” and “second” may explicitly or implicitly include at least one of the features. In the description of the present invention, “plurality” means at least two, for example, two, three, etc. unless expressly specified otherwise.

[0064] For making the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below in conjunction with the drawings and embodiments. It should be understood that the particular embodiments described herein are merely used for explaining the present invention, and are not used for limiting the present invention.

[0065] The present invention is further described now with reference to the drawings of the description.

[0066] Embodiments of the present invention relate to explanation of a method for whole-process numerical simulation and hazard forecast of a mountain disaster. As shown in FIG. 1, the method specifically includes:

[0067] 1. High Space-Time Rainfall Forecast of a Mountain Area

[0068] A precipitation forecast result is comprehensively determined on the basis of computation results of two modes: a global mode (which belongs to a regional rainfall forecast mode having forecast accuracy of 0.125°) of a European center and a mesoscale mode (which belongs to a kilometer-level rainfall dynamic forecast mode) of a southwest area. Firstly, a correction area is determined on the basis of a strong precipitation falling area of the European center, and moreover, a mode computation result of the southwest area is corrected in combination with two pieces of mode forecast data in a delimited precipitation area, so as to improve accuracy. In order to meet the requirement of meteorological data fining required for driving a mountain disaster model, spatial interpolation continue being carried out on a forecast result on the basis of correction, so as to improve space resolution of the data. In order to give consideration to timeliness and objectivity at the same time, a bilinear interpolation method (linear relations are established in longitude and latitude directions respectively, and physical quantity sizes on certain longitudes and latitudes are converted according to the relations) is used for interpolating forecast data having space resolution of 9 KM to 0.01°*0.01° (about 1 km), and then an empirical relation between a physical quantity and an elevation is established according to terrain, so as to describe an influence of the terrain on precipitation.

[0069] The interpolation method uses the bilinear interpolation method, which is also called as bilinear interpolation. The bilinear interpolation is linear interpolation expansion of an interpolation function having two variables, the core idea of which is that one-time linear interpolation is carried out respectively in two directions. An attribute value (FIG. 2) at a point P=(x,y) is computed, linear interpolation is carried out in an x direction to obtain

[00014]$f\left(R_1\right)\approx \frac{x_2-x}{x_2-x_1}f\left(Q_{11}\right)+\frac{x-x_1}{x_2-x_1}f\left(Q_{21}\right),\mathrm{where}$$R_1=\left(x,y_1\right),\mathrm{and}$$f\left(R_2\right)\approx \frac{x_2-x}{x_2-x_1}f\left(Q_{12}\right)+\frac{x-x_1}{x_2-x_1}f\left(Q_{22}\right),\mathrm{where}$$R_2=\left(x,y_2\right),\mathrm{and}$

[0070] then linear interpolation is carried out in a y direction to obtain

[00015]$f\left(P\right)\approx \frac{2_{}-}{2_{}-1_{}}f\left(R_1\right)+\frac{-1_{}}{2_{}-1_{}}f\left(R_2\right)$

[0071] so as to obtain a desired result f(x,):

[00016]$f\left(x,\right)\approx \frac{f\left(Q_{11}\right)}{\left(x_2-x_1\right)\left(2_{}-1_{}\right)}\left(x_2-x\right)\left(2_{}-\right)+\frac{f\left(Q_{21}\right)}{\left(x_2-x_1\right)\left(2_{}-1_{}\right)}\left(x-x_1\right)\left(2_{}-\right)+\frac{f\left(Q_{12}\right)}{\left(x_2-x_1\right)\left(2_{}-1_{}\right)}\left(x_2-x\right)\left(-1_{}\right)+\frac{f\left(Q_{22}\right)}{\left(x_2-x_1\right)\left(2_{}-1_{}\right)}\left(x-x_1\right)\left(-1_{}\right)$

[0072] With reference to FIG. 2, a point P is obtained by data points Q11, Q12, Q22, Q23 and R1 and R2 to be interpolated.

[0073] In the equation, represents a corresponding attribute value at a certain point, and Q.sub.11, Q.sub.12, Q.sub.21 and Q.sub.22 represent point positions at (x.sub.1, y.sub.1), (x.sub.1, y.sub.2), (x.sub.2, y.sub.1) and (x.sub.2, y.sub.2).

[0074] 2. Hydrodynamic Process and Numerical Simulation

(1) Hydrodynamic Process Model

[0075] As shown in FIG. 3, a hydrodynamic module technical route diagram is shown. Depth integral simplification is carried out on the basis of a Naviers-Stokes equation and a convection term in a momentum equation is ignored to obtain a diffusion wave model, and a small watershed hydrodynamic process is quantitatively computed. An Aston vegetation rainfall interception model and a Green-Ampt slope saturation infiltration model are introduced on the basis of the diffusion wave model, and a whole physical event from rainfall to vegetation interception and slope infiltration, and then to runoff generation and motion is simulated, where a governing equation of the model is:

[00017]$\frac{\partial h}{\partial t}+\frac{\partial \mathrm{hu}}{\partial x}+\frac{\partial \mathrm{hv}}{\partial y}=R-V-I$$\frac{\partial \left(1/2{\mathrm{gh}}^2\right)}{\partial x}=-\mathrm{gh}\frac{\partial z_b}{\partial x}+S_{\mathrm{fx}}$$\frac{\partial \left(1/2{\mathrm{gh}}^2\right)}{\partial y}=-\mathrm{gh}\frac{\partial z_b}{\partial y}+S_{\mathrm{fy}}$

[0076] where t is time, x and y are horizontal coordinates, h represents a fluid depth, u and v represent speed components of a fluid depth-averaged speed in the x direction and the y direction respectively, R is reduced rainfall intensity after vegetation interception, I is an infiltration rate, V is vegetation interception, g is gravitational acceleration, Z.sub.b is a substrate elevation, S.sub.fx and S.sub.fy represent substrate friction terms in the x direction and the y direction respectively, and

[0077] S.sub.fx and S.sub.fy are expressed by Manning models as:

[00018]$S_{\mathrm{fx}}=g\frac{n^{2}u\sqrt{u^2+v^2}}{h^{1/3}}$$S_{\mathrm{fy}}=g\frac{n^{2}v\sqrt{u^2+v^2}}{h^{1/3}}$

[0078] where n is a Manning coefficient, h represents the fluid depth, u and v represent the speed components of the fluid depth-averaged speed in the x direction and the y direction respectively, and g is the gravitational acceleration.

[0079] The Aston vegetation rainfall interception model used by vegetation interception V is expressed as:

[00019]$V=S_{\max }\left(1-e^{-k\frac{P_c}{S_{\max }}}\right)$

[0080] where S.sub.max is maximum interception of vegetation, and according to computation of different types of vegetation, P.sub.c is accumulated rainfall, and k is a parameter related to the canopy density of the vegetation, and is expressed as:

k=1−e.sup.−(Co*LAI)

[0081] where Co is the canopy density of the vegetation, and is a leaf area index.

[0082] The infiltration rate I is expressed by the slope saturation infiltration model as:

[00020]$I=\frac{\mathrm{df}}{\mathrm{dt}}=k_s\left[\left({\psi }_f+h\right)\frac{{\theta }_s-{\theta }_i}{f}+1\right]$

[0083] where t is time, k.sub.s is a saturated hydraulic coefficient, ψ.sub.f is a matrix suction water head at a front end of a wetting front, θ.sub.s is saturated moisture content of soil, is initial moisture content of the soil, and f is an accumulated infiltration depth.

(2) Model Solving

[0084] For solving of a two-dimensional diffuse wave equation, in order to simulate a disaster physical process and ensure accuracy, stability and high efficiency of the model and a solution format, numerical solution is carried out by means of a first-order windward difference scheme, and fine-grained parallelization is carried out on a program by utilizing a graphics processing unit (GPU) (graphics card) computation device. There are two main architecture types (NVIDIA and AMD) of the graphics card currently used for high performance computation on the market. Generally, in the same level, the NVIDIA architecture features higher main performance and is more widely applied in various industries earlier, but in recent years, the AMD architecture is also gradually stepped into a high-performance computation mainstream channel, and has strong development momentum. A GPU-like computation card (DCU) is developed on the basis of the AMD architecture in the Sea Light series in China. A system platform has already completed GPU device program schemes of the two mainstream architectures, and different computation devices all have excellent compatibility and computational efficiency for different computation platforms.

(3) Actual Computation Case for Longxi River

[0085] A case research area is a located in a Longhuan river watershed in a northwest mountain area of the Dujiangyan City in Sichuan. The Longxi river is used as a first-level branch of a minjiang river watershed, starts from the Longchi hummock located at the north of the Longxi river, converges into a minjiang river water system from the Duanmu garden located at the south of the Longhua river, and has an overall length of 18.22 km, a watershed area of 79 km.sup.2, and an overall water flow direction from north to south. The Longxi river has branches of Zhucao valley, Lengjin valley, Modao valley, Jianping valley, etc., and numerous valleys. In administrative division, the Longxi river belongs to the Longchi town, and is 17 km away from the Dujiangyan City, the town is bounded by Zipingpu town to the East, Wenchuan county to the west, Zipingpu reservoir to the South and Hongkou township to the north, and the total population is more than 3000. There is Duwen Expressway in the south of Longxi River watershed. Leave from the expressway and take a Longchi tourist route to the Longchi town, or go from an urban area, and transportation is convenient.

[0086] According to rainfall records collected from rainfall data of 9 rainfall stations in a research area from the early morning of June 24 to 23 p.m. of Jun. 26, 2018, and rainfall surface data are prepared by an inverse distance difference method. The rainfall stations and monitoring section positions are shown in FIGS. 4-6 below.

[0087] Due to long duration of rainfall, a total of 9 hours are selected from 20:00 on the 22nd to 5:00 on the 23rd. A simulated water depth result is shown in FIG. 7. By comparing with a measured flow process curve, a result shows that the hydrodynamic model of the system platform may reflect an overall trend of flood, and description of a mountain torrent process basically conforms to a measured result.

[0088] FIGS. 8 and 9 show a vulnerability grading map of a Longxi river watershed and a map of a risk grading result of a mountain disaster in the Longxi river watershed.

[0089] 3. Motion Model and Numerical Simulation of a Mountain Torrent and Debris Flow Disaster

[0090] Due to the advantage of high computation efficiency, surface gravity flow such as flood, high-sand-content water flow and debris flow is usually solved by using a depth averaged shallow water wave equation, including a mass and momentum conservation equations in a motion process of the surface gravity flow. When an erosion or deposition process needs to be considered, an additional equation is needed to reflect a solid phase concentration and an equation for substrate erosion. A depth-averaged continuous medium equation for facies averaging a rainfall-induced debris flow is expressed as:

[00021]$\frac{\partial h}{\partial t}+\frac{\partial \mathrm{hu}}{\partial x}+\frac{\partial \mathrm{hv}}{\partial y}=R-I+\frac{E}{1-p}$$\frac{\partial \mathrm{hu}}{\partial t}+\frac{\partial \left({\mathrm{hu}}^2+0.5{\mathrm{gh}}^2\right)}{\partial x}+\frac{\partial \mathrm{huv}}{\partial y}=-\mathrm{gh}\frac{\partial z_b}{\partial x}-\frac{{\tau }_{\mathrm{fx}}}{\rho }-\frac{\left({\rho }_s-{\rho }_f\right)g{h}^2}{2\rho }\frac{\partial c}{\partial x}+\frac{\left({\rho }_b-\rho \right)\mathrm{uE}}{\rho \left(1-p\right)}$$\frac{\partial \mathrm{hv}}{\partial t}+\frac{\partial \mathrm{huv}}{\partial x}+\frac{\partial \left({\mathrm{hv}}^2+0.5{\mathrm{gh}}^2\right)}{\partial y}=-\mathrm{gh}\frac{\partial z_b}{\partial y}-\frac{{\tau }_{\mathrm{fy}}}{\rho }-\frac{\left({\rho }_s-{\rho }_f\right)g{h}^2}{2\rho }\frac{\partial c}{\partial y}+\frac{\left({\rho }_b-\rho \right)\mathrm{vE}}{\rho \left(1-p\right)}$$\frac{\partial \mathrm{hc}}{\partial t}+\frac{\partial \mathrm{huc}}{\partial x}+\frac{\partial hvc}{\partial y}=E$$\frac{\partial z_b}{\partial t}=-\frac{E}{1-p}$

[0091] where t represents time, x and y are horizontal coordinates, h represents a fluid depth, u and v are speed components of a fluid depth-averaged speed in an x direction and a y direction, respectively, c is a depth-averaged solid phase concentration, g is gravitational acceleration, R is reduced rainfall intensity after vegetation interception, I is an infiltration rate, Z.sub.b is a substrate elevation, ρ is the density of a solid-liquid mixture, and ρ=cρ.sub.s+(1−c)ρ.sub.f, ρ.sub.s and ρ.sub.f are concentrations of a solid phase and a liquid phase respectively, ρ.sub.b is the density of a saturated substrate, ρ.sub.b=(1−p)ρ.sub.s+pρ.sub.f, p is the porosity of a substrate material, τ.sub.fx and τ.sub.fy are substrate resistances in the x direction and the y direction respectively, and E is an erosion rate of the substrate.

[0092] An equation of the erosion rate may be expressed as:

[00022]$E=\frac{{\tau }_f-{\tau }_r}{\rho \sqrt{u^2+v^2}}$

[0093] where shear stress τ.sub.f of the solid-liquid mixture with the substrate may be expressed as a combination τf=cτ.sub.fs÷(1−c)τ.sub.fw of solid phase shear stress and liquid phase shear stress, c is a depth-averaged solid phase concentration, the solid phase shear stress is τ.sub.js=(ρ.sub.s−ρ.sub.f) gh tan ϕ.sub.bed the liquid phase shear stress is τ.sub.fw=ρ.sub.fgn.sup.2(u.sup.2+v.sup.2)h.sup.−1/3, ϕ.sub.bed is a Coulomb friction angle of the substrate, n is a Manning coefficient, ρ.sub.s and ρ.sub.f are concentrations of the solid phase and the liquid phase respectively, h represents the fluid depth, u and v are speed components of the fluid depth-averaged speed along the x direction and y direction respectively, g is gravitational acceleration, the shear stress resistance of a substrate erosion layer may be expressed τ.sub.r=c.sub.o+(1−λ) ρgh tan ϕ.sub.bin, c.sub.o and ϕ.sub.bin are cohesion and an internal friction angle of the substrate material respectively, λ is a pore water pressure coefficient of the substrate, ρ is the density of the solid-liquid mixture, h represents the fluid depth, and g is gravitational acceleration.

[0094] 4. Risk Analysis and Hazard Forecast of a Small Watershed Disaster

(1) Comprehensive Hazard Assessment for the Disaster

[0095] The destructive power of mountain disasters depends on different disaster-causing modes of the disasters, such that hazard degrees of the disasters are different. Considering composite hazard characteristics of impact and burying of mountain disasters, a comprehensive disaster-causing risk degree of the disasters is determined, and a computation model is as follows:

[00023]$H=H_e+H_d$$H_d=\frac{N_{i,j}\Delta V}{A}$$H_e=A.Math.\underset{t>0}{\max }\left[\left(u^2+v^2\right)h\rho \right]$

[0096] where H is the comprehensive disaster-causing risk degree of the disaster, H.sub.e is a hazard caused by impact damage, and is expressed by maximum kinetic energy of a disaster-causing body; H.sup.d is a hazard caused by burying and is represented by a maximum burying depth of the disaster-causing body; N.sub.i,j is equal to the number of particles in a control grid centered on (i,j) and ΔV is the volume of the particle; A is a grid area; and u and v are speeds in the x direction and the y direction respectively, h is the fluid depth, and ρ is the density of the solid-liquid mixture.

(2) Vulnerability Analysis

[0097] The vulnerability analysis of disaster-bearing bodies is mainly used for computing vulnerability of disaster-bearing bodies. In general, the vulnerability of the disaster-bearing bodies is related to values of the disaster-bearing bodies and vulnerability indexes of the disaster-bearing bodies. Vulnerability of the i.sup.th type of disaster-bearing bodies may be described as:

V.sub.i=V(u).sub.i×C.sub.i

[0098] where V.sub.i is vulnerability of the i.sup.th type of disaster-bearing bodies, V(u).sub.t is a comprehensive value of the i.sup.th type of disaster-bearing bodies, and C.sub.i vulnerability indexes of the i.sup.th type of disaster-bearing bodies.

[0099] The quantification of the comprehensive value V(u) of the disaster-bearing bodies depends on an average unit price D.sub.e of the disaster-bearing bodies and an actual area A.sub.e affected by the disaster, and may be expressed as:

V(u)=D.sub.e×A.sub.e

[0100] A value of the vulnerability index C of the disaster-bearing bodies is determined complexly due to the influence of multiple uncertainty factors. Different disaster-bearing bodies are different in positions and damage degrees relative to the debris flow; and the damage degrees of different disaster-bearing bodies are also different when the disaster-bearing bodies are impacted or buried by the same-scale debris flow. The vulnerability index of the disaster-bearing bodies is expressed by a ratio of the fluid depth h to the geometric height H.sup.c of the disaster-bearing body and is expressed as:

[00024]$C=\frac{h}{H_c}$

[0101] where h is a fluid depth, and H.sub.c is a geometric height of the disaster-bearing body, such as the height of a bridge deck, a house, etc. When

[00025]$\frac{h}{H_c}\ge 1,$

it is indicated that the disaster-bearing bodies have been completely buried by the debris flow, and the C is valued as 1.

(3) Risk Assessment and Grading

[0102] Disaster risk evaluation is a comprehensive analysis process about hazard of disaster occurrence and possible harm influence of disaster occurrence, and the risk degree is quantitative expression of the disaster risk. When the disaster occurs, whether the disaster-bearing bodies are exposed to the disaster-causing range of the disaster has uncertainty. On the basis of a numerical simulation result, the risk of the disaster is determined as a comprehensive function of a comprehensive hazard of the disaster, vulnerability of a disaster-bearing body and exposure of the disaster-bearing body, the disaster risk degree is the product of the hazard degree, the vulnerability degree and the exposure degree, a computation equation which is as follows:

[00026]$R=f\left(H,V_d,E_d\right)=\underset{i=1}{\overset{n}{.Math.}}{H}_i\times V_i\times E_i$

[0103] where R is the risk degree of the disaster and is expressed by a certain numerical value ranging from 0 (no risk) to 1 (high risk); H is the comprehensive disaster-causing risk degree of the disaster, V.sub.d is the vulnerability degree of the disaster-bearing body, E.sub.d is the exposure degree of the disaster-bearing body, and H.sub.i is the comprehensive disaster-causing risk degree of the i.sup.th type of disasters, and is expressed by a certain numerical value ranging from 0 (no hazard) to 1 (high hazard); V.sub.i is vulnerability of the i.sup.th type of disaster-bearing bodies and is determined by the damage degree and value or number of the disaster-bearing body; and E.sub.i is the exposure degree of the i.sup.th type of disaster-bearing bodies and is expressed by a certain numerical value ranging from 0 (no exposure) to 1 (complete exposure).

[0104] The above-mentioned embodiments are merely intended for describing the technical solutions of the present invention rather than limiting the present invention. Although the present invention is described in detail with reference to the above-mentioned embodiments, those of ordinary skill in the art should understand that they may still make modifications to the technical solutions described in the embodiments or equivalent substitutions to some or all of the technical features of the technical solutions. These modifications or substitutions do not enable the corresponding technical solutions to depart from the scope of the technical solutions in the embodiments of the present invention, and should all fall with the scope of the claims and the description of the present invention.