CN114491768A - Method for calculating uniform mode production flow of variable production flow layer - Google Patents

Method for calculating uniform mode production flow of variable production flow layer Download PDF

Info

Publication number
CN114491768A
CN114491768A CN202210142924.XA CN202210142924A CN114491768A CN 114491768 A CN114491768 A CN 114491768A CN 202210142924 A CN202210142924 A CN 202210142924A CN 114491768 A CN114491768 A CN 114491768A
Authority
CN
China
Prior art keywords
soil
formula
runoff
intensity
layer
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.)
Granted
Application number
CN202210142924.XA
Other languages
Chinese (zh)
Other versions
CN114491768B (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

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • Y—GENERAL 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
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Geometry (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Algebra (AREA)
  • Data Mining & Analysis (AREA)
  • Structural Engineering (AREA)
  • Civil Engineering (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Architecture (AREA)
  • Operations Research (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Investigation Of Foundation Soil And Reinforcement Of Foundation Soil By Compacting Or Drainage (AREA)

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 in a vertical direction; calculating the evapotranspiration capacity of the computing unit; carrying out runoff yield calculation on the vegetation retention 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 fine description and the rapid forecast of the production convergence mechanism of the distributed hydrological model in the sudden flood forecast in the hilly area are difficult to be considered.

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 uniform mode runoff yield of a variable runoff yield layer.
Background
Most hilly and rocky areas in China are areas where flood burst and easily occur, and due to the comprehensive effects of meteorological and hydrological conditions and special terrains, mountain flood disasters are frequent, so that the mountain flood disasters are difficult points 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 on 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 overall, and some bottlenecks exist in the forecasting of sudden flood in a hilly area. Specifically, the distributed hydrological model not only reflects the physical mechanism of production convergence and 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 for research on 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 raw 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 EpCalculating;
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 2pBased on the water quantity balance principle, an equation is established for the vegetation interception layer of the calculation unit, the interception effect of the calculation unit is represented, and the penetration rain intensity P is calculatedvAnd actual evapotranspiration rate E of vegetation-retaining layerv；
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 principlepAnd penetration rain strength P calculated in step 3vAnd actual evapotranspiration rate EvCalculating the total soil water content M of the field IgEarth-rock interface water level NwtEquivalent constant rain intensity RgSubsurface infiltration strength IgIntensity of interflow qgSaturated surface runoff
Figure BDA0003507162620000021
Intensity of evaporation EgAnd the water flow strength 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 of the step 4 and the water quantity balance principle, the calculation list is subjected toEstablishing an equation in the element soil layer area II, representing the formation and development process of the fluctuation runoff yield layer between soil layers, carrying out runoff yield calculation, and calculating the evapotranspiration capacity E according to the movement wave infiltration principle and the step 2pAnd step 3, calculating the penetration rain strength PvAnd actual evapotranspiration rate EvAnd the subsurface infiltration strength I of the field I calculated in the step 5gWater level N at soil-rock interfacewtCalculating the total soil water content M of the field IIuWater level N between soil layersstEquivalent constant rain intensity RuSubsurface infiltration strength IuWetting front depth NfIntensity of interflow quSaturated surface runoff
Figure BDA0003507162620000023
And evaporation intensity Eu；
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 3vAnd 5, calculating the subsurface infiltration strength I of the field IgAnd saturated surface runoff
Figure BDA0003507162620000024
And the subsurface infiltration strength I of the domain II calculated in step 6uAnd 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 calculatedsrThe replenishment intensity q of the basement rock weathered layer to the soil layerwrAnd groundwater runoff strength qw。
Further, in step 2, the formula of the evapotranspiration capacity of the computing unit is as follows:
Ep＝Ke·Eobs (1)
in the formula (1), KeThe actual measurement evapotranspiration conversion coefficient is obtained; eobsThe evaporation rate was actually observed for the evaporation dish.
Further, in step 3, the water balance equation of the vegetation retaining layer of the unit is calculated as:
Figure BDA0003507162620000031
in the formula (2), MvThe vegetation retention; t is time; p is rainfall intensity; evThe actual evapotranspiration rate of the vegetation-retaining layer; pvPenetration rain intensity; evAnd PvThe corresponding equation is:
Ev＝min(Mv/dt,Ep) (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-weathering 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), IgSubsurface permeability strength for Domain I;
Figure BDA0003507162620000035
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; q. q.swrThe 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.ssr(ii) infiltration makeup strength for domain I to the weathered basement rock formation;
Figure BDA0003507162620000041
saturated surface runoff for domain I formation; egEvaporation intensity for domain I; mgThe 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 the formula (6), NsdIs the soil-weathered bedrock interface depth (soil thickness); n is a radical ofwtThe depth of water level of soil-rock interface is 0-Nwt≤Nsd；θsThe saturated water content of the soil is obtained; r isgIs the equivalent constant rain intensity of domain I; n is the vertical depth of the earth surface (positive when the orientation is downward); theta (R)gAnd 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), K0nThe surface saturation hydraulic conductivity in the n direction; ε ═ 2+3 λ)/λ, where,lambda is the soil pore distribution index; thetarThe 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;
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), it can be deduced
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 notgWhether or not to exceed
Figure BDA0003507162620000048
) The classification into two cases is as follows:
not forming soil-rock interface change runoff layer, i.e.
Figure BDA0003507162620000049
The calculation formula is as follows:
the formula of the water level of the soil-rock interface is as follows:
Nwt＝Nsd (9)
the equivalent constant rain intensity formula is:
Figure BDA0003507162620000051
the subsurface infiltration strength formula is:
Ig＝Rg cosα (11)
the interflow intensity formula is as follows:
qg＝NsdRg(ar-1)sin αL (12)
in the formula (12), arIs the anisotropy ratio; alpha is the gradient of the earth surface; l is the interflow width;
the saturated surface runoff formula is as follows:
Figure BDA0003507162620000052
the evaporation intensity formula is:
Figure BDA0003507162620000053
in the formula (14), the compound represented by the formula (I),
Figure BDA0003507162620000054
is MgBy taking R in the formula (7)gInto equation (6)
Is calculated to obtaingIs empirically taken as 0.01mm/h, EuIs the evaporation intensity of domain II;
the formula of the water flow strength in the n direction at the bottom of the soil layer is as follows:
Figure BDA0003507162620000055
② a soil-rock interface transition runoff layer is formed, namely
Figure BDA0003507162620000056
The calculation formula is as follows:
at the moment, the soil-rock interfaceDepth of water level NwtThe water content of the soil reaches saturation, namely theta (R)g,Nwt)＝θsFrom equations (6) and (7), the soil moisture content formula for domain I can be derived as:
Figure BDA0003507162620000057
formula (16) is NwtAt M, ingAnd 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), c1＝c4(Mg-θsNsd)+c3；c2＝c3c4；c3＝θs-θr；c4-f/epsilon; lambert w function is the inverse of the f (w) w · exp (w) function; the exp function represents an exponential function with a natural constant e as the base; mgThe total soil moisture content of the region I is calculated by the formula (6);
the equivalent constant rain intensity formula is:
Rg＝K0nexp(-fNwt) (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), MuTotal soil moisture in field II; t is time; i isuSubsurface 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 Eu；
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 Nf. Equivalent constant rain intensity R for domain II when the wetting front reaches a certain depthuSaturation hydraulic conductivity above this depth at which interbedded formation of runoff occurs, 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)fWhether or not to exceed
Figure BDA0003507162620000073
) The classification into two cases is as follows:
no 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 soil layers is as follows:
Nst＝Nf (28)
the subsurface infiltration strength formula is:
Iu＝(Ru-Rg)cosα (29)
the interflow intensity formula is as follows:
qu＝Nf(Ru-Rg)(ar-1)sin αL (30)
the saturated surface runoff formula is as follows:
Figure BDA0003507162620000077
the evaporation intensity formula is:
Eu＝min[max(Ep-Ev,0),Mu/dt] (32)
② a mobile runoff layer between soil layers is formed, namely
Figure BDA0003507162620000078
The calculation formula is as follows:
the equivalent constant rain intensity formula is:
Ru＝K0n exp(-fNst) (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 content of the first and second substances,
Figure BDA0003507162620000083
c2＝c3c4；c3＝θs-θr；c4＝-f/ε；
the subsurface infiltration strength formula is:
Figure BDA0003507162620000084
the interflow intensity formula is as follows:
Figure BDA0003507162620000085
the saturated surface runoff strength formula is as follows:
Figure BDA0003507162620000086
the evaporation intensity formula is:
Eu＝min[max(Ep-Ev,0),Mu/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:
qh＝max(Pv-Iu-Ig,0) (40)
Figure BDA0003507162620000087
in the formula (40), qhThe strength of the surface runoff is super-seepage; pvPenetration rain intensity; i isuThe infiltration strength of the soil horizon II; i isgThe infiltration strength of the soil horizon I; in the formula (41), qdSaturated surface runoff strength;
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), MwThe maximum storage capacity of the nonlinear reservoir is
Figure BDA0003507162620000094
Respectively, the ground into which the upstream meshes flowTotal runoff inflow intensity and subsurface runoff intensity flowing to a downstream grid; q. q.ssrThe infiltration strength of the soil layer; q. q.swrThe replenishing strength of 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, the effluent flow strength of the subsurface runoff can be calculated by using a nonlinear motion wave equation considering the degree of fullness, and the formula is as follows:
Figure BDA0003507162620000095
in formula (43), kwIs a weathering bedrock formation outflow parameter; b is a shape parameter; calculating the subsurface runoff intensity of the upstream calculation unit and the underground runoff intensity of the upstream calculation unit, which is converged into the calculation unit, by adopting a formula (43), wherein 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 foundation stratumsrFrom the residual water storage of the weathered basement and the bottom N of the soil layersdIntensity of water flow in n direction
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, the new knowledge of the landslide hydrological experiment is reflected, four runoff components of saturated ground runoff, super-seepage ground runoff, interflow and subsurface runoff can be described at the same time by means of fine depiction of the interflow forming and developing process, the unified description of the full-super-seepage runoff and the conversion mechanism of the full-super-seepage runoff is realized, the artificial assumption of full accumulation or super-seepage of the runoff yield mode is avoided, the problem that the fine description and the quick prediction of the runoff yield mechanism of a distributed hydrological model in the sudden flood forecast of a hilly area are difficult to consider is solved, and the engineering significance is high.
Drawings
FIG. 1 is a flow chart of a method of the present invention for variable runoff layer unified mode runoff computing;
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 humid area, taking a square calculation unit with the area size of 10m by 10m, measuring the soil thickness through a 70-mm-caliber gasoline power drill to determine the depth of a soil-bedrock interface, and dividing the calculation unit into a vegetation retention layer, a soil layer and an weathered bedrock layer according to the vertical direction according to the top end of a canopy (the upper boundary of the calculation unit), a fresh bedrock interface (the lower boundary of the calculation 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 processobsThe evapotranspiration capacity E of the computing unit under sufficient water supplypThe calculation is carried out according to the formula:
Ep＝Ke·Eobs (1)
in the formula (1), KeTaking 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 PvAnd actual evapotranspiration rate E of vegetation-retaining layervThe formula is as follows:
Figure BDA0003507162620000111
in the formula (2), MvThe vegetation retention; t is time; p is rainfall intensity; evThe actual evapotranspiration rate of the vegetation-retaining layer; pvPenetration rain intensity; evAnd PvThe corresponding equation is:
Ev＝min(Mv/dt,Ep) (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 2 mm;
step 4, according to the vertical layering result in the step 1, dynamically dividing a calculation domain of the soil water content theta-surface vertical depth n of the calculation unit soil layer into a domain I and a domain II (shown in a 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 a 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 principlepAnd penetration rain strength P calculated in step 3vAnd actual evapotranspiration rate EvCalculating the total soil moisture content M of the field IgEarth-rock interface water level NwtEquivalent constant rain intensity RgSubsurface infiltration strength IgIntensity of interflow qgSaturated surface runoff
Figure BDA0003507162620000114
Intensity of evaporation EgAnd the strength of the water flow 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), IgSubsurface permeability strength for domain I;
Figure BDA0003507162620000121
are respectively domainsI, total inflow intensity of the lateral interflow of the upstream grid and the intensity of the lateral interflow of the grid area I; a is the area of the calculation unit and is 100m2；qwrThe 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.ssr(ii) infiltration makeup strength for domain I to the weathered basement rock formation;
Figure BDA0003507162620000122
saturated surface runoff for domain I formation; egEvaporation intensity for domain I; mgThe 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 the formula (6), NsdThe soil-weathered bedrock interface depth (soil thickness), determined by the electric drill measurement of step 1, was 1.2 m; n is a radical ofwtThe depth of water level of soil-rock interface is 0-Nwt≤Nsd；θsThe saturated water content of the soil is determined by a drying method and is 0.45; rgIs the equivalent constant rain intensity of domain I; n is the vertical depth of the earth surface (positive when the orientation is downward); theta (R)gAnd 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), K0nThe conductivity of the saturated water power of the earth surface in the n direction is determined by a double-ring measurement method and is 70 mm/h; epsilon is (2+3 lambda)/lambda, wherein lambda is a soil pore distribution index, and the value of epsilon is 4.5; thetarThe 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;
the water content of the soil at the soil-weathered bedrock interface is gradually increased along with the accumulation of the soil water in the region I when the soil thereinAt the critical state of saturation of soil water content (i.e. at
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:
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 notgWhether 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:
Nwt＝Nsd (9)
the equivalent constant rain intensity formula is:
Figure BDA0003507162620000132
the subsurface infiltration strength formula is:
Ig＝Rg cosα (11)
the interflow intensity formula is as follows:
qg＝NsdRg(ar-1)sin αL (12)
in the formula (12), arTaking 20 according to the soil characteristics in the calculation unit according to experience for the anisotropic ratio; alpha is the gradient of the earth surface, and the gradient is 20 degrees and is obtained by carrying out 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 MgBy taking R in the formula (7)gIs calculated by substituting the minimum value of (3) into the formula (6), RgThe minimum value of (A) is empirically taken to be 0.01 mm/h; euIs 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) a soil-rock interface transition runoff layer is formed, namely
Figure BDA0003507162620000137
The calculation formula is as follows:
at the moment, the depth N of the water level of the soil-rock interfacewtThe water content of the soil reaches saturation, namely theta (R)g,Nwt)＝θsFrom equations (6) and (7), the soil moisture content formula for domain I can be derived as:
Figure BDA0003507162620000141
formula (16) is NwtAt M, ingAnd 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), c1＝c4(Mg-θsNsd)+c3；c2＝c3c4；c3＝θs-θr；c4-f/epsilon; lambert w function is the inverse of the function f (w) ═ w · exp (w); mgThe total soil moisture content of the region I is calculated by the formula (6);
the equivalent constant rain intensity formula is:
Rg＝K0n exp(-fNwt) (18)
the subsurface infiltration strength formula is:
Figure BDA0003507162620000143
the interflow intensity formula is as follows:
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 direction of n at the bottom of the soil layer is as follows:
Figure BDA0003507162620000147
step 6, establishing an equation for the soil horizon II of the calculation unit according to the water balance principle, representing the formation and development process of the fluctuation runoff yield layer between soil layers and carrying out runoff yield calculation, and calculating the evapotranspiration capacity E according to the movement wave infiltration principle and the step 2pAnd step 3, calculating the penetration rain strength PvAnd actual evapotranspiration rate EvAnd the subsurface infiltration strength I of the field I calculated in the step 5gWater level N at soil-rock interfacewtCalculating the total soil water content M of the field IIuWater level N between soil layersstEquivalent constant rain intensity RuSubsurface infiltration strength IuWetting front depth NfIntensity of interflow quSaturated surface runoff
Figure BDA0003507162620000151
And evaporation intensity EuAnd calculating a water balance equation of the unit soil horizon II as follows:
Figure BDA0003507162620000152
in formula (24), MuTotal soil moisture in field II; t is time; i isuSubsurface 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 100m2；
Figure BDA0003507162620000154
Saturated surface runoff; evaporation intensity of Domain II is recorded as Eu；
When the penetration rain intensity meets the infiltration supply of the region IThe rest part begins to enter the domain II through infiltration to form a wetting front which continuously develops downwards, and the depth of the wetting front is recorded as Nf. When the wetting front reaches a certain depth, domain II is equivalent to constant rain intensity RuSaturation 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 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 the variation among soil layers (wetting front depth N)fWhether 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 soil layers is as follows:
Nst＝Nf (28)
the subsurface infiltration strength formula is:
Iu＝(Ru-Rg)cosα (29)
the interflow intensity formula is as follows:
qu＝Nf(Ru-Rg)(ar-1)sin αL (30)
the saturated surface runoff formula is as follows:
Figure BDA0003507162620000161
the evaporation intensity formula is:
Eu＝min[max(Ep-Ev,0),Mu/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:
Ru＝K0n exp(-fNst) (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 content of the first and second substances,
Figure BDA0003507162620000165
c2＝c3c4；c3＝θs-θr；c4＝-f/ε；
the subsurface infiltration strength formula is:
Figure BDA0003507162620000166
the interflow intensity formula is as follows:
Figure BDA0003507162620000171
the saturated surface runoff strength formula is as follows:
Figure BDA0003507162620000172
the evaporation intensity formula is:
Eu＝min[max(Ep-Ev,0),Mu/dt] (39)
step 7, according to the vertical layering result of the step 1 and the penetration rain strength P calculated in the step 3vAnd 5, calculating the subsurface infiltration strength I of the field IgAnd saturated surface runoff
Figure BDA0003507162620000173
And the subsurface infiltration strength I of the domain II calculated in step 6uAnd saturated 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:
qh＝max(Pv-Iu-Ig,0) (40)
Figure BDA0003507162620000175
in the formula (40), qhThe strength of the surface runoff is super-seepage; pvPenetration rain intensity; i isuThe infiltration strength of the soil horizon II; i isgThe infiltration strength of the soil horizon I;
in the formula (41), the compound represented by the formula,qdsaturated surface runoff strength;
Figure BDA0003507162620000176
the saturated surface runoff strength formed for the soil horizon II;
Figure BDA0003507162620000177
(ii) the saturated surface runoff intensity formed for soil horizon I;
step 8, calculating the water flow strength of the soil layer bottom in the n direction according to the vertical layering result in the step 1 and the water flow strength of the soil layer bottom 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 strength q of the soil layer is calculatedsrThe replenishment intensity q of the basement rock weathered layer to the soil layerwrAnd groundwater runoff strength qwThe water balance equation is as follows:
Figure BDA0003507162620000179
in formula (42), MwThe 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 basal 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.ssrThe infiltration strength of the soil layer; q. q.swrThe replenishing strength of 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, the effluent flow strength of the subsurface runoff can be calculated by using a nonlinear motion wave equation considering the degree of fullness, and the formula is as follows:
Figure BDA0003507162620000182
in formula (43), kwThe effluence parameter of the weathering bedrock stratum is 0.01 m/h; b is a shape parameter, and the value is 6; calculating the subsurface runoff intensity of the upstream calculation unit and the underground runoff intensity of the upstream calculation unit, which is converged into the calculation unit, by adopting a formula (43), wherein 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 stratumsrFrom the residual water storage of the weathered basement and the bottom N of the soil layersdIntensity 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 supplies the soil layer, and the formula is as follows:
Figure BDA0003507162620000188
a graph of the calculated unit wetting front depth, the water level depth between soil layers and the runoff yield strength 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, calculating unit vertical layering: 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 ispCalculating;
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 2pBased on the water quantity balance principle, an equation is established for the vegetation interception layer of the calculation unit, the interception effect of the calculation unit is represented, and the penetration rain intensity P is calculatedvAnd actual evapotranspiration rate E of vegetation-retaining layerv；
And calculating a water balance equation of the unit vegetation interception layer as follows:
Figure FDA0003507162610000011
in the formula, MvThe vegetation retention; t is time; p is rainfall intensity; evThe actual evapotranspiration rate of the vegetation-retaining layer; pvPenetration 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 a calculation unit soil layer 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 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 calculation unit to represent the formation and development process of the soil-rock interface change runoff formation and development process and carry out runoff formation calculation,
and calculating a water balance equation of the unit soil horizon I as follows:
Figure FDA0003507162610000012
in the formula IgSubsurface permeability strength for Domain I;
Figure FDA0003507162610000013
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; q. q.swrThe 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.ssr(ii) infiltration makeup strength for domain I to the weathered basement rock formation;
Figure FDA0003507162610000021
saturated surface runoff formed for domain I; egEvaporation intensity for domain I; mgTotal water content of soil of domain I;
calculating evapotranspiration capacity E according to the movement wave infiltration principle and the step 2pAnd penetration rain strength P calculated in step 3vAnd actual evapotranspiration rate EvCalculating the total soil water content M of the field IgWater level N of soil-rock interface of region IwtEquivalent constant rain intensity R for domain IgSubsurface infiltration strength of Domain IgIntensity of interflow q of region IgSaturated surface runoff of domain I
Figure FDA0003507162610000022
Evaporation intensity E of Domain IgWater flow intensity in n direction at bottom of soil layer of Heyu I
Figure FDA0003507162610000023
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 FDA0003507162610000024
in the formula, MuTotal soil moisture in field II; t is time; i isuSubsurface permeability strength for domain II;
Figure FDA0003507162610000025
Figure FDA0003507162610000026
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 FDA0003507162610000027
saturated surface runoff for domain II; evaporation intensity of Domain II is recorded as Eu；
Calculating evapotranspiration capacity E according to the movement wave infiltration principle and the step 2pAnd step 3, calculating the penetration rain strength PvAnd actual evapotranspiration rate EvAnd the subsurface infiltration strength I of the field I calculated in the step 5gWater level N at soil-rock interfacewtCalculating the total soil water content M of the field IIuWater level N between soil layers of region IIstEquivalent constant rain intensity R for domain IIuSubsurface infiltration Strength I of Domain IIuRegion II wetting front depth NfIntensity of interflow q of region IIuSaturated surface runoff of domain II
Figure FDA0003507162610000028
And evaporation intensity E of Domain IIu；
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 3vAnd 5, calculating the subsurface infiltration strength I of the field IgAnd saturated surface runoff
Figure FDA0003507162610000029
And the subsurface infiltration strength I of the domain II calculated in step 6uAnd saturated surface runoff
Figure FDA00035071626100000210
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 FDA0003507162610000031
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 calculatedsrThe replenishment intensity q of the basement rock weathered layer to the soil layerwrAnd groundwater runoff strength qw(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 FDA0003507162610000032
in the formula, MwThe maximum storage capacity of the nonlinear reservoir is
Figure FDA0003507162610000033
Figure FDA0003507162610000034
The total inflow intensity of the underground runoff flowing into the upstream grid and the intensity of the underground 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:
Ep＝Ke·Eobs (1)
in the formula, KeThe actual measurement evapotranspiration conversion coefficient is obtained; eobsThe evaporation rate was actually observed for the evaporation dish.
3. The method of claim 1, wherein the method comprises: in said step 5, MgObtained by integrating the soil moisture profile of the region I along the direction n, and the formula is as follows:
Figure FDA0003507162610000035
in the formula, NsdIs the soil-weathered bedrock interface depth; n is a radical ofwtThe depth of water level of soil-rock interface is 0-Nwt≤Nsd；θsThe saturated water content of the soil is obtained; rgIs the equivalent constant rain intensity of domain I; n is the vertical depth of the earth surface; theta (R)gN) unsaturated soil moisture profile distribution for domain I;
Figure FDA0003507162610000036
in the formula, K0nThe surface saturation hydraulic conductivity in the n direction; epsilon is (2+3 lambda)/lambda, wherein lambda is the soil pore distribution index; thetarThe 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 FDA0003507162610000037
Figure FDA0003507162610000041
4. The method of claim 3, wherein the method comprises: 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 FDA0003507162610000042
The formula of the water level of the soil-rock interface is as follows:
Nwt＝Nsd (9)
the equivalent constant rain intensity formula is:
Figure FDA0003507162610000043
the subsurface infiltration strength formula is:
Ig＝Rgcosα (11)
the interflow intensity formula is as follows:
qg＝NsdRg(ar-1)sinαL (12)
in the formula, arIs the anisotropy ratio; alpha is the gradient of the earth surface; l is the interflow width;
the saturated surface runoff formula is as follows:
Figure FDA0003507162610000044
the evaporation intensity formula is:
Figure FDA0003507162610000045
in the formula (I), the compound is shown in the specification,
Figure FDA0003507162610000046
is MgThe minimum value of (d); euIs 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 FDA0003507162610000047
5. the method of claim 3, wherein the method comprises: 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 FDA0003507162610000051
The soil moisture content formula for domain I is:
Figure FDA0003507162610000052
the formula of the water level of the soil-rock interface is as follows:
Figure FDA0003507162610000053
in the formula, c1＝c4(Mg-θsNsd)+c3；c2＝c3c4；c3＝θs-θr；c4-f/epsilon; lambert w function is the inverse of the f (w) w · exp (w) function; mgThe total soil moisture content of the region I is calculated by the formula (6);
the equivalent constant rain intensity formula is:
Rg＝K0nexp(-fNwt) (18)
the subsurface infiltration strength formula is:
Figure FDA0003507162610000054
the interflow intensity formula is as follows:
Figure FDA0003507162610000055
the saturated surface runoff strength formula is as follows:
Figure FDA0003507162610000056
the evaporation intensity formula is:
Figure FDA0003507162610000057
the formula of the water flow strength in the direction of n at the bottom of the soil layer is as follows:
Figure FDA0003507162610000058
Figure FDA0003507162610000067
6. the method of claim 1, wherein the method comprises: in step 6, when the wetting front reaches a certain depth, the equivalent constant rain intensity R of the domain IIuSaturated hydraulic conductivity above the depth at which formation of alternate runoff formation occurs between layers of soil, the depth being measuredIs a critical depth
Figure FDA0003507162610000066
The calculation formula is as follows:
Figure FDA0003507162610000061
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 soil interbed fluctuation is formed or not.
7. The method of claim 6, wherein the method comprises: without formation of a mobile runoff layer between soil layers, i.e.
Figure FDA0003507162610000062
The equivalent constant rain intensity formula is:
Figure FDA0003507162610000063
the wetting front motion formula is:
Figure FDA0003507162610000064
the formula of the depth of the water level between soil layers is as follows:
Nst＝Nf (28)
the subsurface infiltration strength formula is:
Iu＝(Ru-Rg)cosα (29)
the interflow intensity formula is as follows:
qu＝Nf(Ru-Rg)(ar-1)sinαL (30)
the saturated surface runoff formula is as follows:
Figure FDA0003507162610000065
the evaporation intensity formula is:
Eu＝min[max(Ep-Ev,0),Mu/dt] (32)。
8. the method of claim 6, wherein the method comprises: has formed a soil interbedded mobile runoff layer, i.e.
Figure FDA0003507162610000071
The equivalent constant rain intensity formula is:
Ru＝K0nexp(-fNst) (33)
the wetting front motion formula is:
Figure FDA0003507162610000072
the formula of the depth of the water level between soil layers is as follows:
Figure FDA0003507162610000073
in the formula (35), the reaction mixture is,
Figure FDA0003507162610000074
c2＝c3c4；c3＝θs-θr；c4＝-f/ε；
the subsurface infiltration strength formula is:
Figure FDA0003507162610000075
the interflow intensity formula is as follows:
Figure FDA0003507162610000076
the saturated surface runoff strength formula is as follows:
Figure FDA0003507162610000077
the evaporation intensity formula is:
Eu＝min[max(Ep-Ev,0),Mu/dt] (39)。
9. the method of claim 1, wherein the method comprises: 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:
qh＝max(Pv-Iu-Ig,0) (40)
Figure FDA0003507162610000081
in the formula, qhThe strength of the surface runoff is super-seepage; q. q.sdSaturated surface runoff strength.
10. The method of claim 1, wherein the method comprises: in the step 8, 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 accumulation, wherein the formula is as follows:
Figure FDA0003507162610000082
in the formula, kwIs a weathering bedrock formation outflow parameter; 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 FDA0003507162610000083
The underground runoff outflow rate calculated by the calculation unit according to the formula (43) is
Figure FDA0003507162610000084
Infiltration strength q of soil layer to weathered foundation stratumsrFrom the residual water storage of the weathered basement and the bottom N of the soil layersdIntensity of directional water flow
Figure FDA0003507162610000085
Determining, the formula is:
Figure FDA0003507162610000086
when the weathered basement is filled in a non-linear reservoir, i.e.
Figure FDA0003507162610000087
The excess water supplies the soil layer, and the formula is as follows:
Figure FDA0003507162610000088
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 true CN114491768A (en) 2022-05-13
CN114491768B 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)

Cited By (1)

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

Citations (7)

* 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
CN106874605A (en) * 2017-02-22 2017-06-20 中国水利水电科学研究院 A kind of gneiss soil Mountainous Area hillside yardstick hydrologic process analogy method
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

Patent Citations (7)

* 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
CN106874605A (en) * 2017-02-22 2017-06-20 中国水利水电科学研究院 A kind of gneiss soil Mountainous Area hillside yardstick hydrologic process analogy method
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 (2)

* Cited by examiner, † Cited by third party
Title
李彬权等: "基于变动饱和带的产汇流模型及其参数确定方法", 《水科学进展》 *
梁忠民等: "基于统计理论的产流模型", 《水科学进展》 *

Cited By (1)

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

Also Published As

Publication number Publication date
CN114491768B (en) 2022-11-18

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
CN106951612A (en) Dynamic water storage capacity Runoff calculation method in freeze-thawing process of soil
CN109492299A (en) The water resource simulation method coupled based on SWMM with MODFLOW
Schuh et al. Soil moisture redistribution and its effect on inter-annual active layer temperature and thickness variations in a dry loess terrace in Adventdalen, Svalbard
De Louw et al. Rainwater lens dynamics and mixing between infiltrating rainwater and upward saline groundwater seepage beneath a tile-drained agricultural field
CN111563619A (en) Rainfall threshold analysis method for causing watershed landslide risk
CN103306236B (en) Method for constructing underground reservoir in ancient gully of ancient underground river channel
CN114239904A (en) Underground water management method and device
CN114491768B (en) Method for calculating current of variable current layer in unified mode
CN111537698A (en) Rock-soil mass slope surface water movement simulation device with root stone structure and experiment method
CN112257286B (en) Variable-source runoff yield mode simulation method for permafrost region temperature dominance
GUO et al. Rainfall infiltration analysis and infiltration model of slope based on in-situ tests
CN105137042A (en) Method and device for determining construction position of water storage project on Karst slope
CN110286027A (en) Consider that the riverbank soil body of red building root system influence washes away the quantization method of parameter
CN111982780B (en) Method for measuring and calculating relation between urban new district planning underlying surface and old underlying surface
Bahremand et al. Application of WetSpa model for assessing land use impacts on floods in the Margecany–Hornad watershed, Slovakia
CN111999228A (en) Urban new area infiltration measuring and calculating method
Lan et al. A two-dimensional numerical model coupled with multiple hillslope hydrodynamic processes and its application to subsurface flow simulation
Momii et al. Observations and modelling of seawater intrusion for a small limestone island aquifer
Arumí et al. Effect of drought on groundwater in a Chilean irrigated valley
Andréassian et al. Catalogue of the models used in MOPEX 2004/2005
Langsholt Water balance modelling in lateritic terrain
CN213022771U (en) Urban new area infiltration measuring and calculating device
CN108229096A (en) A kind of humid region soil layering Runoff calculation method

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