CN114491768B - Method for calculating current of variable current layer in unified mode - Google Patents

Method for calculating current of variable current layer in unified mode Download PDF

Info

Publication number
CN114491768B
CN114491768B CN202210142924.XA CN202210142924A CN114491768B CN 114491768 B CN114491768 B CN 114491768B CN 202210142924 A CN202210142924 A CN 202210142924A CN 114491768 B CN114491768 B CN 114491768B
Authority
CN
China
Prior art keywords
formula
soil
runoff
layer
intensity
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210142924.XA
Other languages
Chinese (zh)
Other versions
CN114491768A (en
Inventor
赵建飞
梁忠民
段雅楠
李彬权
仇知雨
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN202210142924.XA priority Critical patent/CN114491768B/en
Publication of CN114491768A publication Critical patent/CN114491768A/en
Application granted granted Critical
Publication of CN114491768B publication Critical patent/CN114491768B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/10Geometric CAD
    • G06F30/13Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/28Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00Details relating to the application field
    • G06F2113/08Fluids
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Abstract

The invention discloses a method for calculating the variable production flow layer unified mode production flow, which comprises the following steps: layering the computing units vertically; calculating the evapotranspiration capacity of the computing unit; carrying out runoff calculation on the vegetation interception layer; dynamically dividing a soil layer theta-n calculation domain into a domain I and a domain II according to a steady-state soil moisture profile, respectively establishing an equation for the domain I and the domain II, representing the formation and development process of a soil-rock interface variable runoff generating layer and a soil interlayer variable runoff generating layer, and calculating the runoff generating; the method is used for calculating the runoff yield of the weathered basal rock stratum, and realizes the unified description of the accumulation-super-seepage runoff yield and the conversion mechanism thereof by finely depicting the formation and development process of the interflow runoff. The method solves the problem that the production convergence mechanism of the distributed hydrological model in the sudden flood forecast of the hilly area is difficult to be finely described and quickly forecasted.

Description

Method for calculating uniform mode production flow of variable production flow layer
Technical Field
The invention belongs to the technical field of hydrological models, and particularly relates to a method for calculating a current in a uniform mode of a variable current layer.
Background
Most hilly areas in China are areas with sudden and easy flood occurrence, and due to the comprehensive action of meteorological and hydrological conditions and special terrains, mountain flood disasters are frequent, so that the mountainous areas are difficult and weak links for flood prevention work in China. The hydrological model is a main tool for flood forecasting, and can be generally divided into a black box model, a conceptual model and a distributed hydrological model. Research and application show that a data-oriented black box/conceptual model is difficult to meet the requirement of sudden flood forecasting in a hilly area, and a distributed hydrological model becomes a future trend because the distributed hydrological model can realize flood forecasting at any position/section, and the trend is a necessary direction for the development of the hydrological model. Therefore, the research of the distributed hydrological model facing to the hilly area has important significance.
Since 1969 when the FH69 blueprint is proposed, a plurality of distributed hydrological models which are based on different theories and have different definition and complexity in the hydrological process are appeared, but the progress is gradually slowed down on the whole, and some bottlenecks exist in the sudden flood forecast in a hilly area. Specifically, the distributed hydrological model not only reflects the physical mechanism of production convergence but also is sufficiently refined, but also meets the requirement of rapid forecasting of sudden flood in the hilly area, so that a model structure, particularly a production flow calculation method, needs to be broken through urgently, and the distributed hydrological model is a key problem in research of the hilly area distributed hydrological model.
Disclosure of Invention
The invention aims to overcome the technical defects and provide a method for calculating the yield of a varied yield layer in a unified mode, so as to solve the problem that the detailed description and the rapid forecast of a yield convergence mechanism of a distributed hydrological model in the sudden flood forecast in a hilly area are difficult to be considered at the same time.
In order to solve the technical problems, the technical scheme adopted by the invention is as follows:
a method for calculating the current of a variable current layer in a unified mode comprises the following steps after a calculation unit is determined through original data acquisition:
step 1, calculating unit vertical layering: for a specific computing unit, vertically dividing the computing unit into a vegetation retention layer, a soil layer and a weathered bed according to the top of a canopy (an upper boundary of the computing unit), a fresh bed rock interface (a lower boundary of the computing unit), the earth surface and the soil-bed rock interface;
step 2, calculating the evaporation capacity of the unit: according to the vertical layering result in the step 1, assuming that the evapotranspiration occurs in a vegetation retention layer and a soil layer, the evapotranspiration under the condition of fully supplying water to the computing unit is the evapotranspiration capacity E p Calculating;
step 3, calculating the runoff yield of the vegetation interception layer: according toStep 1, vertical layering result, rainfall intensity P and evapotranspiration capacity E calculated in step 2 p Based on the principle of water balance, an equation is established for the vegetation interception layer of the computing unit, the interception function of the computing unit is represented, and the penetration rain strength P is computed v And actual evapotranspiration rate E of vegetation-retaining layer v
Step 4, dividing a soil layer theta-n calculation domain: according to the vertical layering result in the step 1, dynamically dividing a calculation domain of the soil water content theta-earth surface vertical depth n of a calculation unit soil layer into a domain I and a domain II by adopting a stable soil water profile;
step 5, calculating the runoff yield of the soil horizon I: according to the division result in the step 4 and the water quantity balance principle, an equation is established for the soil horizon I of the calculation unit, the formation and development processes of the soil-rock interface change runoff formation are represented, runoff calculation is carried out, and the evapotranspiration capacity E calculated in the step 2 is calculated according to the movement wave infiltration principle p And penetration rain strength P calculated in step 3 v And actual evapotranspiration rate E v Calculating the total soil water content M of the field I g Earth-rock interface water level N wt Equivalent constant rain intensity R g Subsurface infiltration strength I g Intensity of interflow q g Saturated surface runoff
Figure BDA0003507162620000021
Intensity of evaporation E g And the strength of the water flow in the direction of n at the bottom of the soil layer
Figure BDA0003507162620000022
Step 6, calculating the runoff yield of the soil horizon II: according to the division result in the step 4 and the water quantity balance principle, an equation is established for the soil horizon II of the calculation unit, the formation and development process of the fluctuation runoff formation among the soil horizon is represented, the runoff formation calculation is carried out, and the evapotranspiration capacity E calculated in the step 2 is calculated according to the movement wave infiltration principle p And 3, calculating the penetration rain strength P v And actual evapotranspiration rate E v And the subsurface infiltration strength I of the field I calculated in the step 5 g Water level N at soil-rock interface wt Calculating the total soil water content M of the field II u Soil interlayer waterBit N st Equivalent constant rain intensity R u Subsurface infiltration strength I u Wetting front depth N f Intensity of interflow q u Saturated surface runoff
Figure BDA0003507162620000023
And evaporation intensity E u
Step 7, calculating the soil stratum runoff yield: according to the vertical layering result of the step 1 and the penetration rain strength P calculated in the step 3 v And 5, calculating the subsurface infiltration strength I of the field I g And saturated surface runoff
Figure BDA0003507162620000024
And the calculated subsurface infiltration intensity I of the region II in step 6 u And saturated surface runoff
Figure BDA0003507162620000025
Solving the super-seepage surface runoff and the saturated surface runoff of the computing unit based on a runoff formation principle;
step 8, calculating the weathered basal rock stratum production flow: according to the vertical layering result in the step 1 and the water flow intensity in the direction of n at the bottom of the soil layer calculated in the step 5
Figure BDA0003507162620000026
Based on the water balance principle, an equation is established for the weathering bedrock layer of the calculation unit, the storage and discharge effects of the weathering bedrock in the calculation unit are represented, and the infiltration strength q of the soil layer is calculated sr The replenishment intensity q of the basement rock weathered layer to the soil layer wr And groundwater runoff strength q w
Further, in step 2, the formula of the evapotranspiration capacity of the computing unit is as follows:
E p =K e ·E obs (1)
in the formula (1), K e The actual measurement evapotranspiration conversion coefficient is obtained; e obs The evaporation rate was actually observed for the evaporation dish.
Further, in the step 3, the water balance equation of the unit vegetation entrapment layer is calculated as follows:
Figure BDA0003507162620000031
in the formula (2), M v The vegetation retention; t is time; p is rainfall intensity; e v The actual evapotranspiration rate of the vegetation-retaining layer; p v Penetration rain intensity; e v And P v The corresponding equation is:
E v =min(M v /dt,E p ) (3)
Figure BDA0003507162620000032
in the formula (4), the reaction mixture is,
Figure BDA0003507162620000033
the maximum cut-off for vegetation.
Further, in the step 4, a calculation domain for calculating the soil water content θ of the soil layer of the unit soil layer and the vertical depth n of the earth surface is dynamically divided into two parts, namely a domain I and a domain II according to the steady-state soil water profile, wherein the upper boundary of the domain I is the earth surface, the lower boundary of the domain I is a soil-weathered bedrock interface, the upper boundary of the domain II is the earth surface, and the lower boundary of the domain II is a soil-rock interface water level. At the initial moment, the water level of the soil-rock interface is superposed with the soil-weathered bedrock interface; along with the accumulation of moisture, a 'soil-rock interface change runoff layer' is gradually formed on the soil-weathering bedrock interface, and the water level of the soil-rock interface begins to rise; with the continuous rainfall, the water level of the soil-rock interface is continuously lifted until the water level reaches the earth surface, the area II disappears, the whole soil layer is classified as an area I and is saturated;
further, in step 5, the water balance equation of the unit soil horizon I is calculated as:
Figure BDA0003507162620000034
in the formula (5), I g Subsurface permeability strength for domain I;
Figure BDA0003507162620000035
respectively obtaining total inflow intensity of the lateral subsurface flow of the grid at the upstream of the area I and the lateral subsurface flow intensity of the area I; a is the area of the calculation unit; q. q.s wr The replenishment strength of the weathered basal rock layer to the region I after the weathered basal rock layer reaches the maximum storage capacity; q. q.s sr The infiltration replenishment strength of the domain I to the weathered basement formation;
Figure BDA0003507162620000041
saturated surface runoff for domain I formation; e g Evaporation intensity for domain I; m is a group of g The total soil moisture content of the region I is obtained by integrating the soil moisture profile of the region I along the direction n, and the formula is as follows:
Figure BDA0003507162620000042
in formula (6), N sd Is the soil-weathered bedrock interface depth (soil thickness); n is a radical of wt The depth of water level of soil-rock interface is 0-N wt ≤N sd ;θ s The saturated water content of the soil is obtained; r g Is the equivalent constant rain intensity of domain I; n is the vertical depth of the earth's surface (positive under orientation); theta (R) g And n) unsaturated soil moisture profile distribution of the domain I, according to the theory of the infiltration of the motion wave, comprising:
Figure BDA0003507162620000043
in the formula (7), K 0n The surface saturation hydraulic conductivity in the n direction; epsilon = (2 +3 lambda)/lambda, wherein lambda is a soil pore distribution index; theta.theta. r The residual water content of the soil; e is a natural constant; f is the attenuation coefficient of the saturation hydraulic conductivity along with the depth;
as the soil moisture in the region I accumulates, the soil moisture content at the soil-weathered bedrock interface gradually increases, when the soil moisture content reaches a critical state of saturation (i.e., the soil moisture content reaches a critical state of saturation)
Figure BDA0003507162620000044
) Forming a soil-rock interface change runoff layer, wherein the corresponding soil water content total critical value is
Figure BDA0003507162620000045
From the equations (6) and (7) can be derived
Figure BDA0003507162620000046
The formula of (1) is:
Figure BDA0003507162620000047
wherein the exp function represents an exponential function with a natural constant e as a base;
when the runoff yield calculation is carried out according to the water balance equation of the formula (5), the runoff yield layer (the total soil water content M) is changed according to whether the soil-rock interface is formed or not g Whether or not to exceed
Figure BDA0003507162620000048
) The method is divided into two cases, specifically as follows:
(1) without formation of altered stratosphere at the earth-rock interface, i.e.
Figure BDA0003507162620000049
The calculation formula is as follows:
the formula of the water level of the soil-rock interface is as follows:
N wt =N sd (9)
the equivalent constant rain intensity formula is:
Figure BDA0003507162620000051
the subsurface infiltration strength formula is:
I g =R g cosα (11)
the interflow intensity formula is as follows:
q g =N sd R g (a r -1)sin αL (12)
in the formula (12), a r Is the anisotropy ratio; alpha is the gradient of the earth surface; l is the interflow width;
the saturated surface runoff formula is:
Figure BDA0003507162620000052
the evaporation intensity formula is:
Figure BDA0003507162620000053
in the formula (14), the compound represented by the formula (I),
Figure BDA0003507162620000054
is M g By taking R in the formula (7) g Into equation (6)
Is calculated to obtain g Is empirically taken to be 0.01mm/h, E u Is the evaporation intensity of domain II;
the formula of the water flow strength in the direction of n at the bottom of the soil layer is as follows:
Figure BDA0003507162620000055
(2) a soil-rock interface transition runoff layer is formed, namely
Figure BDA0003507162620000056
The calculation formula is as follows:
at the moment, the depth N of the water level of the soil-rock interface wt Reaches saturation of water content of soil, i.e. theta (R) g ,N wt )=θ s From equations (6) and (7), the soil moisture content formula for domain I can be derived as:
Figure BDA0003507162620000057
formula (16) is N wt At M, in g And when other parameters are known, the formula can be reversely solved according to the formula, and then the formula of the water level of the soil-rock interface is as follows:
Figure BDA0003507162620000058
in the formula (17), c 1 =c 4 (M gs N sd )+c 3 ;c 2 =c 3 c 4 ;c 3 =θ sr ;c 4 = f/epsilon; the LambertW function is an inverse function of f (w) = w.exp (w) function; the exp function represents an exponential function with a natural constant e as the base; m g The total soil moisture content of the region I is calculated by the formula (6);
the equivalent constant rain intensity formula is:
R g =K 0n exp(-fN wt ) (18)
the subsurface infiltration strength formula is:
Figure BDA0003507162620000061
the interflow intensity formula is as follows:
Figure BDA0003507162620000062
the saturated surface runoff strength formula is as follows:
Figure BDA0003507162620000063
the evaporation intensity formula is:
Figure BDA0003507162620000064
the formula of the water flow strength in the direction of n at the bottom of the soil layer is as follows:
Figure BDA0003507162620000065
further, in step 6, the water balance equation of the unit soil horizon II is calculated as:
Figure BDA0003507162620000066
in formula (24), M u Total soil moisture in field II; t is time; i is u Subsurface permeability strength for domain II;
Figure BDA0003507162620000067
respectively obtaining the total inflow intensity of the lateral interflow of the upstream grid of the region II and the lateral interflow intensity of the local grid region II; a is the area of the calculation unit;
Figure BDA0003507162620000068
saturated surface runoff; evaporation intensity of Domain II is recorded as E u
When the penetration rain strength meets the infiltration supply of the domain I, the rest part enters the domain II through infiltration to form a wetting front which continuously develops downwards, and the depth of the wetting front is recorded as N f . Equivalent constant rain intensity R for domain II when the wetting front reaches a certain depth u Saturation hydraulic conductivity above this depth at which a variable runoff formation occurs between layers of soil, this depth being referred to as the critical depth
Figure BDA0003507162620000071
The calculation formula is as follows:
Figure BDA0003507162620000072
when the runoff yield calculation is performed according to the equation (24), the runoff yield layer is formed according to the variation among soil layers (wetting front depth N) f Whether or not to exceed
Figure BDA0003507162620000073
) The classification into two cases is as follows:
(1) without formation of a mobile runoff layer between soil layers, i.e.
Figure BDA0003507162620000074
The calculation formula is as follows:
the equivalent constant rain intensity formula is:
Figure BDA0003507162620000075
the wetting front motion formula is:
Figure BDA0003507162620000076
the formula of the depth of the water level between the soil layers is as follows:
N st =N f (28)
the subsurface infiltration strength formula is:
I u =(R u -R g )cosα (29)
the interflow intensity formula is as follows:
q u =N f (R u -R g )(a r -1)sin αL (30)
the saturated surface runoff formula is:
Figure BDA0003507162620000077
the evaporation intensity formula is:
E u =min[max(E p -E v ,0),M u /dt] (32)
(2) has formed a soil interbedded mobile runoff layer, i.e.
Figure BDA0003507162620000078
Formula for calculationThe following:
the equivalent constant rain intensity formula is:
R u =K 0n exp(-fN st ) (33)
the wetting front motion formula is:
Figure BDA0003507162620000081
the formula of the depth of the water level between the soil layers is as follows:
Figure BDA0003507162620000082
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0003507162620000083
c 2 =c 3 c 4 ;c 3 =θ sr ;c 4 =-f/ε;
the subsurface infiltration strength formula is:
Figure BDA0003507162620000084
the interflow intensity formula is:
Figure BDA0003507162620000085
the saturated surface runoff strength formula is as follows:
Figure BDA0003507162620000086
the evaporation intensity formula is:
E u =min[max(E p -E v ,0),M u /dt] (39)
further, in step 7, the formula of the super-seepage surface runoff and the formula of the saturated surface runoff of the soil layer are respectively as follows:
q h =max(P v -I u -I g ,0) (40)
Figure BDA0003507162620000087
in the formula (40), q h The strength of the surface runoff is super-seepage; p v Penetration rain intensity; i is u The infiltration strength of soil horizon II; i is g The infiltration strength of the soil horizon I; in the formula (41), q d Saturated surface runoff intensity;
Figure BDA0003507162620000091
the saturated surface runoff strength formed for the soil horizon II;
Figure BDA0003507162620000092
(ii) the saturated surface runoff intensity formed for soil horizon I;
further, in step 8, a storage and discharge process of the unit weathering basement rock stratum is calculated by adopting a nonlinear reservoir generalization, and a water balance equation is as follows:
Figure BDA0003507162620000093
in formula (42), M w The maximum storage capacity of the nonlinear reservoir is
Figure BDA0003507162620000094
The total inflow intensity of the subsurface runoff flowing into the upstream grid and the intensity of the subsurface runoff flowing into the downstream grid are respectively; q. q.s sr The infiltration strength of the soil layer; q. q.s wr The replenishing strength to the soil layer when the nonlinear reservoir is fully stored; assuming that the gradient of the fresh bedrock interface is consistent with the gradient of the earth surface, calculating the runoff outflow intensity of the underground runoff by using a nonlinear moving wave equation considering the degree of fullness, wherein the formula is as follows:
Figure BDA0003507162620000095
in formula (43), k w Is a weathering bedrock formation outflow parameter; b is a shape parameter; the underground runoff intensity of the upstream computing unit imported into the computing unit is computed by adopting a formula (43), and the sum is
Figure BDA0003507162620000096
The underground runoff outflow intensity calculated by the calculation unit according to the formula (43) is
Figure BDA0003507162620000097
Infiltration strength q of soil layer to weathered basement rock stratum sr From the residual water storage of the weathered basement and the bottom N of the soil layer sd Intensity of water flow in n direction of the water
Figure BDA0003507162620000098
Determining, the formula is:
Figure BDA0003507162620000099
when the weathered basement is filled in a non-linear reservoir, i.e.
Figure BDA00035071626200000910
The excess water supplies the soil layer, and the formula is as follows:
Figure BDA00035071626200000911
the prior art is referred to in the art for techniques not mentioned in the present invention.
The invention achieves the following beneficial effects: according to the method for calculating the runoff yield of the varied runoff yield layer in the unified mode, the description of the interflow is expanded to the soil-rock interface from a single soil layer by introducing the soil-weathered bedrock interface, new knowledge of a landslide hydrological experiment is reflected, four runoff components of saturated surface runoff, super-seepage surface runoff, interflow and underground runoff can be described at the same time by means of fine depiction of the interflow forming and developing process, unified description of the full-accumulation-super-seepage runoff yield and a conversion mechanism of the full-accumulation-super-seepage yield is realized, artificial assumption of full accumulation or super-seepage of a runoff mode is avoided, the problem that fine description and rapid prediction of a runoff yield mechanism of a distributed hydrological model in hilly area sudden flood prediction are difficult to take into account is solved, and the method has strong engineering significance.
Drawings
FIG. 1 is a flow chart of a variable runoff layer unified mode runoff computing method of the present invention;
FIG. 2 is a schematic view of a computing unit being vertically layered;
FIG. 3 is a schematic diagram of the division of the theta-n domain of the soil layer of the calculation unit;
FIG. 4 is a graph of the unit wetting front depth, the water level depth between soil layers and the runoff yield strength over time as calculated from a slope in the southern wetting area.
Detailed Description
The invention is further described below with reference to the accompanying drawings. The following examples are only for illustrating the technical solutions of the present invention more clearly, and the protection scope of the present invention is not limited thereby.
As shown in figure 1 of the drawings, in which,
step 1, for a hillside with good vegetation in a southern wet area, taking a square computing unit with the area size of 10m × 10m, measuring the soil thickness through a 70-mm-caliber gasoline power drill to determine the depth of a soil-bedrock interface, and vertically dividing the computing unit into a vegetation retaining layer, a soil layer and an weathered bedrock layer according to the top end of a canopy (the upper boundary of the computing unit), a fresh bedrock interface (the lower boundary of the computing unit), the earth surface and the soil-bedrock interface (figure 2);
step 2, according to the vertical layering result in the step 1, assuming that evapotranspiration occurs in a vegetation retaining layer and a soil layer, and taking corresponding evaporation pan observation data E in the rainfall process obs The evapotranspiration capacity E of the computing unit under sufficient water supply p The calculation is carried out according to the formula:
E p =K e ·E obs (1)
in the formula (1), K e Taking 0.97 for actually measured evapotranspiration conversion coefficient;
step 3, according to the vertical layering result in the step 1 and the evapotranspiration capacity calculated in the step 2, taking the actually measured rainfall data corresponding to the rainfall process in the step 2 and converting the actually measured rainfall data into rainfall intensity P, establishing a water balance equation for the vegetation retention layer of the calculation unit based on the water balance principle, representing the retention effect of the calculation unit, and calculating the penetration rainfall intensity P v And actual evapotranspiration rate E of vegetation-retaining layer v The formula is as follows:
Figure BDA0003507162620000111
in the formula (2), M v The vegetation retention; t is time; p is rainfall intensity; e v The actual evapotranspiration rate of the vegetation-retaining layer; p is v Penetration rain intensity; e v And P v The corresponding equation is:
E v =min(M v /dt,E p ) (3)
Figure BDA0003507162620000112
in the formula (4), the reaction mixture is,
Figure BDA0003507162620000113
determining the maximum interception amount of vegetation by calculating the vegetation condition of a unit, and taking 2mm;
step 4, according to the vertical layering result in the step 1, dynamically dividing a calculation domain of the soil water content theta-the surface vertical depth n of the calculation unit soil layer into a domain I and a domain II (shown in figure 3) by adopting a stable soil water profile, wherein the upper boundary of the domain I is the surface of the earth, the lower boundary of the domain I is a soil-weathered bedrock interface, the upper boundary of the domain II is the surface of the earth, the lower boundary of the domain II is a soil-rock interface water level, and at the initial moment, the soil-rock interface water level is superposed with the soil-weathered bedrock interface (shown in figure 3-a); as moisture accumulates, a "soil-rock interface transition runoff layer" gradually forms on the soil-weathering bedrock interface, and the soil-rock interface water level begins to rise (fig. 3-b and 3-c); with the continuous rainfall, the water level of the soil-rock interface is lifted continuously until the water level reaches the earth surface, the area II disappears, the whole soil layer is classified as an area I and is saturated (fig. 3-d);
step 5, according to the division result of the step 4 and the water quantity balance principle, an equation is established for the soil horizon I of the calculation unit, the formation and development processes of the soil-rock interface change runoff formation are represented, runoff calculation is carried out, and the evapotranspiration capacity E calculated in the step 2 is calculated according to the movement wave infiltration principle p And penetration rain strength P calculated in step 3 v And actual evapotranspiration rate E v Calculating the total soil moisture content M of the field I g Earth-rock interface water level N wt Equivalent constant rain intensity R g Subsurface infiltration strength I g Intensity of interflow q g Saturated surface runoff
Figure BDA0003507162620000114
Intensity of evaporation E g And the water flow strength in the direction of n at the bottom of the soil layer
Figure BDA0003507162620000115
And calculating a water balance equation of the unit soil horizon I as follows:
Figure BDA0003507162620000116
in the formula (5), I g Subsurface permeability strength for Domain I;
Figure BDA0003507162620000121
respectively obtaining the total inflow intensity of the lateral interflow of the grid at the upstream of the region I and the intensity of the lateral interflow of the grid region I; a is the area of the calculation unit and is 100m 2 ;q wr The replenishment strength of the weathered basal rock layer to the region I after the weathered basal rock layer reaches the maximum storage capacity; q. q of sr The infiltration replenishment strength of the domain I to the weathered basement formation;
Figure BDA0003507162620000122
formed for domain IAnd surface runoff; e g Evaporation intensity for domain I; m is a group of g The total soil moisture content of the region I is obtained by integrating the soil moisture profile of the region I along the direction n, and the formula is as follows:
Figure BDA0003507162620000123
in formula (6), N sd The soil-weathered bedrock interface depth (soil thickness), determined by the electric drill measurement of step 1, was 1.2m; n is a radical of hydrogen wt The depth of water level of soil-rock interface is 0-N wt ≤N sd ;θ s The saturated water content of the soil is determined by a drying method and is 0.45; r is g Is the equivalent constant rain intensity of domain I; n is the vertical depth of the earth's surface (positive in orientation); theta (R) g And n) unsaturated soil moisture profile distribution of the domain I, according to the theory of the infiltration of the motion wave, comprising:
Figure BDA0003507162620000124
in the formula (7), K 0n The conductivity of the saturated water power of the earth surface in the n direction is determined by a double-ring measurement method and is 70mm/h; epsilon = (2 +3 lambda)/lambda, wherein lambda is a soil pore distribution index, and the value of epsilon is 4.5; theta r The residual water content of the soil is determined by a drying method and is 0.093; f is a saturated hydraulic conductivity attenuation coefficient along with depth, and is taken according to the soil profile of the calculation unit according to experience;
as the soil moisture in the region I accumulates, the soil moisture content at the soil-weathered bedrock interface gradually increases, when the soil moisture content reaches a critical state of saturation (i.e., the soil moisture content reaches a critical state of saturation)
Figure BDA0003507162620000125
) Forming a soil-rock interface change runoff layer, wherein the corresponding soil water content total critical value is
Figure BDA0003507162620000126
From the equations (6) and (7), it can be deduced
Figure BDA0003507162620000127
The formula of (1) is as follows:
Figure BDA0003507162620000128
when the runoff yield calculation is carried out according to the water balance equation of the formula (5), the runoff yield layer (the total soil water content M) is changed according to whether the soil-rock interface is formed or not g Whether or not to exceed
Figure BDA0003507162620000129
) The classification into two cases is as follows:
1) Without formation of altered stratosphere at the earth-rock interface, i.e.
Figure BDA0003507162620000131
The calculation formula is as follows:
the formula of the water level of the soil-rock interface is as follows:
N wt =N sd (9)
the equivalent constant rain intensity formula is:
Figure BDA0003507162620000132
the subsurface infiltration strength formula is:
I g =R g cosα (11)
the interflow intensity formula is as follows:
q g =N sd R g (a r -1)sin αL (12)
in the formula (12), a r Taking 20 according to the soil characteristics in the calculation unit according to experience if the ratio is anisotropic; alpha is the gradient of the earth surface, and the gradient is 20 degrees, which is obtained by performing gradient analysis on the DEM; l is the subsurface flow width and is 10m obtained by DEM analysis;
the saturated surface runoff formula is as follows:
Figure BDA0003507162620000133
the evaporation intensity formula is:
Figure BDA0003507162620000134
in the formula (14), the compound represented by the formula (I),
Figure BDA0003507162620000135
is M g By taking R in the formula (7) g Is calculated by substituting the minimum value of (3) into the formula (6), R g The minimum value of (2) is empirically taken as 0.01mm/h; e u Is the evaporation intensity of domain II.
The formula of the water flow strength in the direction of n at the bottom of the soil layer is as follows:
Figure BDA0003507162620000136
2) Has formed a soil-rock interface transition runoff layer, i.e.
Figure BDA0003507162620000137
The calculation formula is as follows:
at the moment, the depth N of the water level of the soil-rock interface wt The water content of the soil reaches saturation, namely theta (R) g ,N wt )=θ s From equations (6) and (7), the soil moisture content formula for domain I can be derived as:
Figure BDA0003507162620000141
formula (16) is N wt At M, in g And when other parameters are known, the formula can be reversely solved according to the formula, and then the formula of the water level of the soil-rock interface is as follows:
Figure BDA0003507162620000142
in the formula (17), c 1 =c 4 (M gs N sd )+c 3 ;c 2 =c 3 c 4 ;c 3 =θ sr ;c 4 = -f/epsilon; the LambertW function is an inverse function of f (w) = w.exp (w) function; m g The total soil moisture content of the region I is calculated by the formula (6);
the equivalent constant rain intensity formula is:
R g =K 0n exp(-fN wt ) (18)
the subsurface infiltration strength formula is:
Figure BDA0003507162620000143
the interflow intensity formula is:
Figure BDA0003507162620000144
the saturated surface runoff strength formula is as follows:
Figure BDA0003507162620000145
the evaporation intensity formula is:
Figure BDA0003507162620000146
the formula of the water flow strength in the n direction at the bottom of the soil layer is as follows:
Figure BDA0003507162620000147
step 6, establishing an equation for the soil horizon II of the computing unit according to the water balance principle, and representing the formation and development process of the variable runoff yield layer between the soil horizonAnd calculating the runoff yield, and calculating the evapotranspiration capacity E according to the movement wave infiltration principle and the step 2 p And step 3, calculating the penetration rain strength P v And actual evapotranspiration rate E v And the subsurface infiltration intensity I of the field I calculated in step 5 g And the water level N of the soil-rock interface wt Calculating the total soil water content M of the field II u Water level N between soil layers st Equivalent constant rain intensity R u Subsurface infiltration strength I u Wetting front depth N f Intensity of interflow q u Saturated surface runoff
Figure BDA0003507162620000151
And evaporation intensity E u And calculating a water balance equation of the unit soil horizon II as follows:
Figure BDA0003507162620000152
in formula (24), M u Total soil moisture in Domain II; t is time; i is u Subsurface permeability strength for domain II;
Figure BDA0003507162620000153
respectively obtaining the total inflow intensity of the lateral interflow of the upstream grid of the region II and the lateral interflow intensity of the local grid region II; a is the area of the calculation unit and is 100m 2
Figure BDA0003507162620000154
Saturated surface runoff; evaporation intensity of Domain II is recorded as E u
When the penetration rain strength meets the infiltration supply of the domain I, the rest part enters the domain II through infiltration to form a wetting front which continuously develops downwards, and the depth of the wetting front is recorded as N f . When the wetting front reaches a certain depth, the equivalent constant rain intensity R of the domain II u Saturation hydraulic conductivity above this depth at which interbedded formation of runoff occurs, this depth being referred to as the critical depth
Figure BDA0003507162620000155
The calculation formula is as follows:
Figure BDA0003507162620000156
when the runoff yield calculation is performed according to the equation (24), the runoff yield layer is formed according to whether the soil interbed fluctuation is generated (wetting front depth N) f Whether or not to exceed
Figure BDA0003507162620000157
) The classification into two cases is as follows:
1) Without formation of a mobile runoff layer between soil layers, i.e.
Figure BDA0003507162620000158
The calculation formula is as follows:
the equivalent constant rain intensity formula is:
Figure BDA0003507162620000159
the wetting front motion formula is:
Figure BDA00035071626200001510
the formula of the depth of the water level between the soil layers is as follows:
N st =N f (28)
the subsurface infiltration strength formula is:
I u =(R u -R g )cosα (29)
the interflow intensity formula is as follows:
q u =N f (R u -R g )(a r -1)sin αL (30)
the saturated surface runoff formula is as follows:
Figure BDA0003507162620000161
the evaporation intensity formula is:
E u =min[max(E p -E v ,0),M u /dt] (32)
2) Has formed a soil interbedded mobile runoff layer, i.e.
Figure BDA0003507162620000162
The calculation formula is as follows:
the equivalent constant rain intensity formula is:
R u =K 0n exp(-fN st ) (33)
the wetting front motion formula is:
Figure BDA0003507162620000163
the formula of the depth of the water level between soil layers is as follows:
Figure BDA0003507162620000164
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0003507162620000165
c 2 =c 3 c 4 ;c 3 =θ sr ;c 4 =-f/ε;
the subsurface infiltration strength formula is:
Figure BDA0003507162620000166
the interflow intensity formula is:
Figure BDA0003507162620000171
the saturated surface runoff strength formula is as follows:
Figure BDA0003507162620000172
the evaporation intensity formula is:
E u =min[max(E p -E v ,0),M u /dt] (39)
step 7, according to the vertical layering result of the step 1 and the penetration rain strength P calculated in the step 3 v And 5, calculating the subsurface infiltration strength I of the field I g And saturation of surface runoff
Figure BDA0003507162620000173
And the calculated subsurface infiltration intensity I of the region II in step 6 u And saturation of surface runoff
Figure BDA0003507162620000174
Based on the runoff formation principle, solving the super-seepage surface runoff and the saturated surface runoff of the calculation unit, wherein the super-seepage surface runoff and the saturated surface runoff of the soil layer are respectively represented by the following formulas:
q h =max(P v -I u -I g ,0) (40)
Figure BDA0003507162620000175
in the formula (40), q h The strength of the surface runoff is super-seepage; p is v Penetration rain intensity; I.C. A u The infiltration strength of soil horizon II; I.C. A g The infiltration strength of the soil horizon I;
in the formula (41), q d Saturated surface runoff strength;
Figure BDA0003507162620000176
the saturated surface runoff intensity formed for the soil horizon II;
Figure BDA0003507162620000177
the saturated surface runoff intensity formed for the soil horizon I;
step 8, calculating the water flow intensity in the n direction at the bottom of the soil layer according to the vertical layering result in the step 1 and the water flow intensity in the n direction at the bottom of the soil layer calculated in the step 5
Figure BDA0003507162620000178
Based on the water balance principle, a nonlinear reservoir equation is established for the weathering bedrock of the calculation unit, the storage and discharge effects of the weathering bedrock in the calculation unit are represented, and the infiltration intensity q of the soil layer is calculated sr And the supplement strength q of the bedrock weathered layer to the soil layer wr And groundwater runoff strength q w The water balance equation is as follows:
Figure BDA0003507162620000179
in formula (42), M w The maximum storage capacity of the nonlinear reservoir is
Figure BDA00035071626200001710
According to the geological condition of the calculation unit and the development condition of the weathered basement rock stratum, taking 0.8m according to experience;
Figure BDA0003507162620000181
the total inflow intensity of the subsurface runoff flowing into the upstream grid and the intensity of the subsurface runoff flowing into the downstream grid are respectively; q. q of sr The infiltration strength of the soil layer; q. q of wr The replenishing strength to the soil layer when the nonlinear reservoir is fully stored; assuming that the gradient of the fresh bedrock interface is consistent with the gradient of the earth surface, calculating the runoff outflow intensity of the underground runoff by using a nonlinear moving wave equation considering the degree of fullness, wherein the formula is as follows:
Figure BDA0003507162620000182
in formula (43), k w The value of the efflorescence bedrock outflow parameter is 0.01m/h; b is a shape parameter, and the value is 6; the underground runoff intensity of the upstream computing unit imported into the computing unit is computed by adopting a formula (43), and the sum is
Figure BDA0003507162620000183
The underground runoff outflow intensity calculated by the calculation unit according to the formula (43) is
Figure BDA0003507162620000184
Infiltration strength q of soil layer to weathered foundation stratum sr From the residual water storage of the weathered basement and the bottom N of the soil layer sd Intensity of water flow in n direction
Figure BDA0003507162620000185
Determining, the formula is:
Figure BDA0003507162620000186
when the weathered basement is filled in a non-linear reservoir, i.e.
Figure BDA0003507162620000187
The excess water is supplied to the soil layer by the formula:
Figure BDA0003507162620000188
a graph of the calculated unit wetting front depth, the water level depth between soil layers and the runoff intensity over time is shown in FIG. 4.
The above description is only a preferred embodiment of the present invention, and it should be noted that, for those skilled in the art, several modifications and variations can be made without departing from the technical principle of the present invention, and these modifications and variations should also be regarded as the protection scope of the present invention.

Claims (10)

1. A method for calculating the current of a variable current layer in a unified mode is characterized by comprising the following steps:
step 1, vertically layering a computing unit: dividing the computing unit into a vegetation interception layer, a soil layer and an weathered bed rock layer according to the top end of the canopy, the fresh bed rock interface, the earth surface and the soil-bed rock interface in the vertical direction;
step 2, calculating the evaporation capacity of the unit: according to the vertical layering result of the step 1, the evapotranspiration is assumed to occur in a vegetation retaining layer and a soil layer, and the evapotranspiration capacity E is p Calculating;
step 3, calculating the runoff yield of the vegetation interception layer: according to the vertical layering result in the step 1, the rainfall intensity P and the evapotranspiration capacity E calculated in the step 2 p Based on the principle of water balance, an equation is established for the vegetation interception layer of the computing unit, the interception function of the computing unit is represented, and the penetration rain strength P is computed v And actual evapotranspiration rate E of vegetation-retaining layer v
And calculating a water balance equation of the unit vegetation interception layer as follows:
Figure FDA0003885928930000011
in the formula, M v The vegetation retention; t is time; p is rainfall intensity; e v The actual evapotranspiration rate of the vegetation-retaining layer; p v Penetration rain intensity;
step 4, dividing a soil layer theta-n calculation domain: according to the vertical layering result in the step 1, dynamically dividing a calculation domain of the soil water content theta-earth surface vertical depth n of the soil layer of the calculation unit into a domain I and a domain II by adopting a stable soil water profile; wherein the upper boundary of the region I is the earth surface, the lower boundary is a soil-weathered bedrock interface, the upper boundary of the region II is also the earth surface, and the lower boundary is a soil-rock interface water level;
step 5, calculating the runoff yield of the soil horizon I: according to the division result in the step 4 and the water balance principle, an equation is established for the soil horizon I of the computing unit to represent the formation and development process of the soil-rock interface fluctuation runoff formation and the runoff formation calculation,
and calculating a water balance equation of the unit soil horizon I as follows:
Figure FDA0003885928930000012
in the formula I g Subsurface permeability strength for Domain I;
Figure FDA0003885928930000013
respectively obtaining total inflow intensity of the lateral subsurface flow of the grid at the upstream of the area I and the lateral subsurface flow intensity of the area I; a is the area of the calculation unit; q. q.s wr The replenishment strength of the weathered basal rock layer to the region I after the weathered basal rock layer reaches the maximum storage capacity; q. q of sr The infiltration replenishment strength of the domain I to the weathered basement formation;
Figure FDA0003885928930000014
saturated surface runoff formed for domain I; e g Evaporation intensity for domain I; m is a group of g Total soil moisture of Domain I; t is time;
calculating evapotranspiration capacity E according to the movement wave infiltration principle and the step 2 p And penetration rain strength P calculated in step 3 v And actual evapotranspiration rate E v Calculating the total soil water content M of the field I g Water level N at soil-rock interface of region I wt Equivalent constant rain intensity R for domain I g Subsurface infiltration Strength I of Domain I g Intensity of interflow q of region I g Saturated surface runoff of domain I
Figure FDA0003885928930000021
Evaporation intensity E of Domain I g Water flow intensity in the n direction at the bottom of soil layer of Heyu I
Figure FDA0003885928930000022
Step 6, calculating the runoff yield of the soil horizon II: according to the division result of the step 4 and the water balance principle, an equation is established for the soil horizon II of the calculation unit to represent the formation and development process of the fluctuation runoff yield layer among the soil horizon and carry out runoff yield calculation,
and calculating a water balance equation of the unit soil horizon II as follows:
Figure FDA0003885928930000023
in the formula, M u Total soil moisture in Domain II; t is time; i is u Subsurface permeability strength for domain II;
Figure FDA0003885928930000024
Figure FDA0003885928930000025
respectively obtaining the total inflow intensity of the lateral interflow of the upstream grid of the region II and the lateral interflow intensity of the local grid region II; a is the area of the calculation unit;
Figure FDA0003885928930000026
saturated surface runoff for domain II; evaporation intensity of Domain II is recorded as E u
Calculating evapotranspiration capacity E according to the movement wave infiltration principle and the step 2 p And 3, calculating the penetration rain strength P v And actual evapotranspiration rate E v And the subsurface infiltration intensity I of the field I calculated in step 5 g Water level N at soil-rock interface wt Calculating the total soil water content M of the field II u Water level N between soil layers of region II st Equivalent constant rain intensity R for domain II u Subsurface infiltration Strength I of Domain II u Region II wetting front depth N f Intensity of interflow q of region II u Saturated surface runoff of region II
Figure FDA0003885928930000027
And evaporation intensity E of field II u
Step 7, calculating the soil stratum runoff yield: according to the vertical layering result of the step 1 and the penetration rain strength P calculated in the step 3 v And 5, calculating the subsurface infiltration strength I of the field I g And saturation of surface runoff
Figure FDA0003885928930000028
And the subsurface infiltration strength I of the domain II calculated in step 6 u And to the saturationSurface runoff
Figure FDA0003885928930000029
Solving the super-seepage surface runoff and the saturated surface runoff of the computing unit based on a runoff formation principle;
step 8, weathered basement rock stratum runoff yield calculation: according to the vertical layering result in the step 1 and the water flow intensity in the direction of n at the bottom of the soil layer calculated in the step 5
Figure FDA00038859289300000210
Based on a water balance principle, an equation is established for the weathering bedrock layer of the calculation unit, the storage and discharge effects of the weathering bedrock in the calculation unit are represented, and the infiltration replenishment strength q of the domain I to the weathering bedrock layer is calculated sr Replenishment strength q of the weathered basement rock formation to the zone I after reaching its maximum storage wr And groundwater runoff strength q w (ii) a The storage and discharge process of the unit weathering basement rock stratum is calculated by adopting nonlinear reservoir generalization, and the water balance equation is as follows:
Figure FDA0003885928930000031
in the formula, M w Is the storage capacity of a nonlinear reservoir, and the maximum storage capacity is
Figure FDA0003885928930000032
Figure FDA0003885928930000033
The total inflow intensity of the subsurface runoff flowing into the upstream grid and the intensity of the subsurface runoff flowing into the downstream grid are respectively.
2. The method of claim 1, wherein the method comprises:
in the step 2, the formula of the evapotranspiration capacity of the computing unit is as follows:
E p =K e ·E obs (1)
in the formula, K e The actual measurement evapotranspiration conversion coefficient is obtained; e obs The evaporation rate was actually observed for the evaporation dish.
3. The method of claim 1 for calculating a uniform pattern of a variable stream layer, wherein:
in said step 5, M g Obtained by integrating the soil moisture profile of the region I along the direction n, and the formula is as follows:
Figure FDA0003885928930000034
in the formula, N sd Is the soil-weathered bedrock interface depth; n is a radical of wt The depth of water level of soil-rock interface is 0-N wt ≤N sd ;θ s The saturated water content of the soil is obtained; r is g Is the equivalent constant rain intensity for domain I; n is the vertical depth of the earth surface; theta (R) g N) unsaturated soil moisture profile distribution for domain I;
Figure FDA0003885928930000035
in the formula, K 0n The surface saturation hydraulic conductivity in the n direction; ε = (2 +3 λ)/λ, where λ is the soil pore distribution index; theta r The residual water content of the soil; e is a natural constant; f is the attenuation coefficient of the saturated hydraulic conductivity along with the depth;
when the water content of the soil reaches a saturated critical state, a soil-rock interface change runoff layer is formed, and the corresponding critical value of the total water content of the soil is
Figure FDA0003885928930000036
Figure FDA0003885928930000037
4. The method of claim 3 for calculating a uniform pattern of a variable stream layer, wherein: carrying out runoff yield calculation according to the water balance equation shown in the formula (5) when no soil-rock interface change runoff yield layer is formed, namely
Figure FDA0003885928930000041
The formula of the water level of the soil-rock interface is as follows:
N wt =N sd (9)
the equivalent constant rain intensity formula is:
Figure FDA0003885928930000042
the subsurface infiltration strength formula is:
I g =R g cosα (11)
the interflow intensity formula is:
q g =N sd R g (a r -1)sinαL (12)
in the formula, a r Is the anisotropy ratio; alpha is the gradient of the earth surface; l is the interflow width;
the saturated surface runoff formula is:
Figure FDA0003885928930000043
the evaporation intensity formula is:
Figure FDA0003885928930000044
in the formula (I), the compound is shown in the specification,
Figure FDA0003885928930000045
is M g Minimum value of (d); e u Is the evaporation intensity of domain II;
the formula of the water flow strength in the direction of n at the bottom of the soil layer is as follows:
Figure FDA0003885928930000046
5. the method of claim 3 for calculating a uniform pattern of a variable stream layer, wherein: carrying out runoff formation calculation according to the water balance equation shown in the formula (5) to form a soil-rock interface change runoff formation, namely
Figure FDA0003885928930000047
The soil moisture content formula for domain I is:
Figure FDA0003885928930000048
the formula of the water level of the soil-rock interface is as follows:
Figure FDA0003885928930000051
in the formula, c 1 =c 4 (M gs N sd )+c 3 ;c 2 =c 3 c 4 ;c 3 =θ sr ;c 4 = -f/epsilon; the LambertW function is an inverse function of f (w) = w.exp (w) function; m g The total soil moisture content of the region I is calculated by the formula (6);
the equivalent constant rain intensity formula is:
R g =K 0n exp(-fN wt ) (18)
the subsurface infiltration strength formula is:
Figure FDA0003885928930000052
in the formula, alpha is the gradient of the earth surface;
the interflow intensity formula is as follows:
Figure FDA0003885928930000053
in the formula, a r Is the anisotropy ratio; l is the subsurface flow width;
the saturated surface runoff strength formula is as follows:
Figure FDA0003885928930000054
the evaporation intensity formula is:
Figure FDA0003885928930000055
in the formula (I), the compound is shown in the specification,
Figure FDA0003885928930000056
is M g The minimum value of (d);
the formula of the water flow strength in the n direction at the bottom of the soil layer is as follows:
Figure FDA0003885928930000057
6. the method of claim 1 for calculating a uniform pattern of a variable stream layer, wherein: in step 6, when the wetting front reaches a certain depth, the equivalent constant rain intensity R of the domain II u Saturation hydraulic conductivity above this depth at which interbedded formation of runoff occurs, this depth being referred to as the critical depth
Figure FDA0003885928930000058
It calculatesThe formula is as follows:
Figure FDA0003885928930000061
in the formula, f is the attenuation coefficient of saturated hydraulic conductivity along with depth; k 0n The surface saturation hydraulic conductivity in the n direction;
when the runoff yield calculation is performed according to the equation (24), the runoff yield layer is divided into two cases according to whether the runoff yield layer is formed between the soil layers.
7. The method of claim 6, wherein the method comprises: without formation of a mobile runoff layer between soil layers, i.e.
Figure FDA0003885928930000062
The equivalent constant rain intensity formula is:
Figure FDA0003885928930000063
in the formula, epsilon = (2 +3 lambda)/lambda, wherein lambda is a soil pore distribution index; theta.theta. s The saturated water content of the soil is obtained; theta r The residual water content of the soil;
the wetting front motion formula is:
Figure FDA0003885928930000064
in the formula, alpha is the gradient of the earth surface;
the formula of the depth of the water level between soil layers is as follows:
N st =N f (28)
the subsurface infiltration strength formula is:
I u =(R u -R g )cosα (29)
the interflow intensity formula is as follows:
q u =N f (R u -R g )(a r -1)sinαL (30)
in the formula, a r Is the anisotropy ratio; l is the subsurface flow width;
the saturated surface runoff formula is as follows:
Figure FDA0003885928930000071
the evaporation intensity formula is:
E u =min[max(E p -E v ,0),M u /dt] (32)。
8. the method of claim 6 for calculating a uniform pattern of a variable stream layer, wherein: has formed a variable runoff layer between soil layers, i.e.
Figure FDA0003885928930000072
The equivalent constant rain intensity formula is:
R u =K 0n exp(-fN st ) (33)
the wetting front motion formula is:
Figure FDA0003885928930000073
in the formula, alpha is the gradient of the earth surface; theta s The saturated water content of the soil is obtained;
the formula of the depth of the water level between soil layers is as follows:
Figure FDA0003885928930000074
in the formula (35), the reaction mixture is,
Figure FDA0003885928930000075
c 2 =c 3 c 4 ;c 3 =θ sr ;c 4 = f/epsilon; the LambertW function is an inverse function of f (w) = w.exp (w) function; theta (R) g N) unsaturated soil moisture profile distribution for domain I; epsilon = (2 +3 lambda)/lambda, wherein lambda is a soil pore distribution index; theta.theta. r The residual water content of the soil;
the subsurface infiltration strength formula is:
Figure FDA0003885928930000076
the interflow intensity formula is as follows:
Figure FDA0003885928930000081
in the formula, a r The anisotropy ratio is adopted, and L is the interflow flow width;
the saturated surface runoff intensity formula is as follows:
Figure FDA0003885928930000082
the evaporation intensity formula is:
E u =min[max(E p -E v ,0),M u /dt] (39)。
9. the method of claim 1 for calculating a uniform pattern of a variable stream layer, wherein: in the step 7, the formula of the super-seepage surface runoff and the formula of the saturated surface runoff of the soil layer are respectively as follows:
q h =max(P v -I u -I g ,0) (40)
Figure FDA0003885928930000083
in the formula, q h The surface runoff intensity is super-seepage; q. q of d Saturated surface runoff strength.
10. The method of claim 1 for calculating a uniform pattern of a variable stream layer, wherein: in the step 8, assuming that the gradient of the fresh bedrock interface is consistent with the gradient of the earth surface, calculating the subsurface runoff strength by using a nonlinear motion wave equation considering the degree of fullness, wherein the formula is as follows:
Figure FDA0003885928930000084
in the formula, k w An effluent parameter for the weathered basement formation; b is a shape parameter; the underground runoff quantity of the upstream computing unit and the underground runoff quantity of the upstream computing unit which are converged into the computing unit are calculated by adopting a formula (43), and the sum is
Figure FDA0003885928930000085
The underground runoff outflow rate calculated by the calculation unit according to the formula (43) is
Figure FDA0003885928930000086
Alpha is the gradient of the earth surface;
infiltration recharge strength q of Domain I to Weathered Foundation sr The residual water storage capacity of the weathered basal rock stratum and the water flow strength of the n direction at the bottom of the soil layer of the region I
Figure FDA0003885928930000087
Determining, the formula is:
Figure FDA0003885928930000091
when the weathered bed is not linearly reservoir full, i.e.
Figure FDA0003885928930000092
The excess water supplies the soil layer, and the formula is as follows:
Figure FDA0003885928930000093
CN202210142924.XA 2022-02-16 2022-02-16 Method for calculating current of variable current layer in unified mode Active CN114491768B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210142924.XA CN114491768B (en) 2022-02-16 2022-02-16 Method for calculating current of variable current layer in unified mode

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210142924.XA CN114491768B (en) 2022-02-16 2022-02-16 Method for calculating current of variable current layer in unified mode

Publications (2)

Publication Number Publication Date
CN114491768A CN114491768A (en) 2022-05-13
CN114491768B true CN114491768B (en) 2022-11-18

Family

ID=81481733

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210142924.XA Active CN114491768B (en) 2022-02-16 2022-02-16 Method for calculating current of variable current layer in unified mode

Country Status (1)

Country Link
CN (1) CN114491768B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115176549B (en) * 2022-09-14 2023-03-28 中国煤炭地质总局勘查研究总院 Ecological restoration method for seasonal frozen soil in plateau alpine region

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102419788A (en) * 2010-12-16 2012-04-18 南京大学 Method for designing distributed-type hydrographical model based on penetration-storage integrated dynamic runoff yield mechanism
CN106599605A (en) * 2017-02-22 2017-04-26 中国水利水电科学研究院 Method for simulating hydrologic process of hillside scale in limestone earth-rock mountain area
CN106951612A (en) * 2017-03-06 2017-07-14 河海大学 Dynamic water storage capacity Runoff calculation method in freeze-thawing process of soil
CN108229096A (en) * 2018-03-13 2018-06-29 河海大学 A kind of humid region soil layering Runoff calculation method
WO2021217776A1 (en) * 2020-04-27 2021-11-04 中山大学 Hydro-meteorological time sequence-based runoff simulation method and system
AU2021106591A4 (en) * 2021-05-21 2021-11-18 Chinese Research Academy Of Environmental Sciences Method And System For Calculating A Dynamic Water Environment Capacity Of A Basin Based On River Lake Water Quality Limit

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106874605B (en) * 2017-02-22 2018-12-28 中国水利水电科学研究院 A kind of gneiss soil Mountainous Area hillside scale hydrologic process analogy method

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102419788A (en) * 2010-12-16 2012-04-18 南京大学 Method for designing distributed-type hydrographical model based on penetration-storage integrated dynamic runoff yield mechanism
CN106599605A (en) * 2017-02-22 2017-04-26 中国水利水电科学研究院 Method for simulating hydrologic process of hillside scale in limestone earth-rock mountain area
CN106951612A (en) * 2017-03-06 2017-07-14 河海大学 Dynamic water storage capacity Runoff calculation method in freeze-thawing process of soil
CN108229096A (en) * 2018-03-13 2018-06-29 河海大学 A kind of humid region soil layering Runoff calculation method
WO2021217776A1 (en) * 2020-04-27 2021-11-04 中山大学 Hydro-meteorological time sequence-based runoff simulation method and system
AU2021106591A4 (en) * 2021-05-21 2021-11-18 Chinese Research Academy Of Environmental Sciences Method And System For Calculating A Dynamic Water Environment Capacity Of A Basin Based On River Lake Water Quality Limit

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于变动饱和带的产汇流模型及其参数确定方法;李彬权等;《水科学进展》;20211220;第208-218页 *

Also Published As

Publication number Publication date
CN114491768A (en) 2022-05-13

Similar Documents

Publication Publication Date Title
Fu et al. Role of epikarst in near‐surface hydrological processes in a soil mantled subtropical dolomite karst slope: Implications of field rainfall simulation experiments
CN108776851B (en) Method for determining early warning threshold value of shallow landslide disaster induced by rainstorm
CN109492299A (en) The water resource simulation method coupled based on SWMM with MODFLOW
CN106951612A (en) Dynamic water storage capacity Runoff calculation method in freeze-thawing process of soil
CN103306236B (en) Method for constructing underground reservoir in ancient gully of ancient underground river channel
De Louw et al. Rainwater lens dynamics and mixing between infiltrating rainwater and upward saline groundwater seepage beneath a tile-drained agricultural field
CN114491768B (en) Method for calculating current of variable current layer in unified mode
CN114239904A (en) Underground water management method and device
CN106640076A (en) Water-retention coal mining method for cooperatively controlling water level and water quantity of unconfined aquifer
CN107916646A (en) A kind of rain-flood resources amount of taping the latent power evaluation method
CN105160121A (en) Finite element control based modeling method for distributed hydrological model
CN112257286B (en) Variable-source runoff yield mode simulation method for permafrost region temperature dominance
CN106202911A (en) A kind of river channel ecology flow recharge of ground water needs water computational methods
CN105137042A (en) Method and device for determining construction position of water storage project on Karst slope
Bahremand et al. Application of WetSpa model for assessing land use impacts on floods in the Margecany–Hornad watershed, Slovakia
CN109085023B (en) Karst region rock-soil interface flow efficient collection method and device
Zhang et al. Numerical simulation of groundwater recharge from south-to-north water division project
Arumí et al. Effect of drought on groundwater in a Chilean irrigated valley
Momii et al. Observations and modelling of seawater intrusion for a small limestone island aquifer
Lan et al. A two-dimensional numerical model coupled with multiple hillslope hydrodynamic processes and its application to subsurface flow simulation
CN111999228A (en) Urban new area infiltration measuring and calculating method
Yuan et al. Insights into the design and development of shanghai coastal reservoirs
Hoxie Digital model of the Arikaree aquifer near Wheatland, southeastern Wyoming
CN114429089B (en) Distributed nonlinear hydrological simulation method for karst area
CN107085632A (en) Underground engineering water buoyancy computational methods under the conditions of weak/impermeable stratum

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant