Efficient water eutrophication evaluation method integrating inorganic and organic substances

Abstract

The invention discloses an efficient water eutrophication evaluation method integrating inorganic and organic substances. The method comprises the following steps: firstly, acquiring DEM topographic data as an input condition of a two-dimensional hydrodynamic force-substance transportation coupling model; assigning values to related parameters in the model by using the discrete DEM topographic data of the rectangular grid; then constructing a substance transport diffusion equation; constructing a two-dimensional hydrodynamic force-substance transportation coupling control equation set; the concentration gradient of each environmental factor in the model is limited; converting a biochemical reaction item Ssch in a source item of the hydrodynamic force-substance transport coupling equation set into each substance state variable conversion balance equation for one-by-one solving and the like, and finally evaluating the eutrophication condition of the water body through a comprehensive index method. According to the method, the occurrence mechanism of the eutrophication of the river and the lake is clarified in detail by comprehensively considering the organic and inorganic forms of nitrogen and phosphorus substances in the water body, so that the change of the eutrophication state of the water body and the improvement degree of the corresponding water quality by implementation of engineering measures are effectively evaluated.

Description

Technical Field

[0001] The invention belongs to the field of water environment and water ecology technology in rivers, lakes and seas, and specifically relates to a high-efficiency water eutrophication assessment method that integrates inorganic and organic substances. Background Art

[0002] Eutrophication of water bodies leads to the massive growth of algae, which reduces the transparency of water bodies, causes black and smelly phenomena, and causes hypoxia, resulting in the death of a large number of fish. For this reason, the simulation and prediction of eutrophication of water bodies has become the focus of attention of many experts and scholars, and is also one of the hot spots in the field of water environment and water ecology. Many scholars are committed to developing water eutrophication models, and single nutrient load models and phytoplankton and nutrient salt related models have emerged. Such as the Vollenweider model, the Droop model, the estimation of phytoplankton primary productivity, and the ecological dynamics model (WASP model). Compared with the previous models, the ecological dynamics model can simulate the eutrophication of water bodies in more detail and accurately. However, this type of model involves many ecosystem state variables and the transformation relationship between them is relatively complex, resulting in low model simulation efficiency. Based on this, an efficient water eutrophication assessment method that integrates inorganic and organic substances is constructed. By comprehensively considering the organic and inorganic forms of nitrogen and phosphorus substances in water bodies, the material, chemical and biological processes inside water bodies are mastered, and the occurrence mechanism of river and lake eutrophication is explained. In this way, the effectiveness of the implementation of engineering measures in changing the eutrophic status of water bodies and the corresponding water quality improvements is evaluated. Summary of the invention

[0003] The purpose of the present invention is to provide an efficient water eutrophication assessment method that integrates inorganic and organic substances. By comprehensively considering the organic and inorganic forms of nitrogen and phosphorus substances in water bodies, the occurrence mechanism of river and lake eutrophication is explained in detail, so as to effectively evaluate the changes in the eutrophication status of water bodies and the corresponding improvement in water quality caused by the implementation of engineering measures.

[0004] The technical solution adopted by the present invention is a high-efficiency water eutrophication assessment method that integrates inorganic and organic substances, which is specifically implemented according to the following steps:

[0005] Step 1: Obtain DEM terrain data as input conditions for the two-dimensional hydrodynamic-material transport coupling model;

[0006] Step 2: Use the rectangular grid discretized DEM terrain data to assign values ​​to the relevant parameters in the model;

[0007] Step 3: Construct the material transport and diffusion equation;

[0008] Step 4: construct a two-dimensional hydrodynamic-material transport coupling control equation set;

[0009] Step 5: Limit the concentration gradient of each environmental factor in the model;

[0010] Step 6: Substitute the biochemical reaction term S in the source term of the hydrodynamic-material transport coupling equation sc h is converted into the state variable transformation equilibrium equation of each substance and solved one by one;

[0011] Step 7: define the bottom slope elevation of the unit grid boundary surface as the maximum value at the boundaries of two adjacent units for reconstruction, and reconstruct the corresponding water depth and water level;

[0012] Step 8: Achieve the spatiotemporal second-order accuracy of the model;

[0013] Step 9: Evaluate the eutrophication of rivers and lakes;

[0014] Step 10: Evaluate the eutrophication status of water bodies through the comprehensive index method.

[0015] The present invention is also characterized in that:

[0016] The material transport and diffusion equation in step 3 is expressed as follows:

[0017]

[0018] In the formula, C s is the average substance concentration along the water depth direction, in kg / m 3 ;D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s; S in source and sink items sc h is the biochemical reaction term, and its unit is kg / m 3 ; C0 is the initial average substance concentration along the water depth direction, in kg / m 3 .

[0019] The vector form of the two-dimensional hydrodynamic-mass transport coupling equations in step 4 is as follows:

[0020]

[0021] in,

[0022]

[0023] Where t is time in seconds; Q is a variable vector, where h is water depth in meters, and q is x and q y They are the single width flow in the x and y directions, in m 2 / s; F, G are flux vectors in the x and y directions respectively; u, v are flow velocities in the x and y directions respectively, in units of m / s; S is the source term vector, including bottom slope source term, friction source term, material diffusion term and material source term; z b is the riverbed bottom elevation, in m; C f is the friction coefficient of the bed surface, calculated as C f =gn 2 / h 1 / 3 , where n is the Manning coefficient; g is the acceleration due to gravity, in m / s 2 ; C s is the average substance concentration along the water depth direction, in kg / m 3 ;D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s; S in source and sink items sc h is the biochemical reaction term, and its unit is kg / m 3 ; C0 is the initial average substance concentration along the water depth direction, in kg / m 3 .

[0024] Step 4 is as follows:

[0025] Step 4.1: In the control unit grid i, the central finite volume method is used to integrate the hydrodynamic-material transport diffusion coupling control equations. Based on the Gaussian divergence theorem, the area integral of the flux term is converted into a line integral in the form of:

[0026]

[0027] Where Ω is the volume of control volume i; Г is the boundary of control volume i; n is the unit vector in the direction of the external normal corresponding to the control volume boundary Г; F(Q)·n is the external normal convection term flux vector on the control volume boundary Г; S is the source term vector, including the bottom slope source term, friction source term, diffusion term and biochemical reaction term,

[0028] The flux F(Q)·n line integral is:

[0029]

[0030] Where Ω is the volume of the i-th control unit grid, in m 3 ; Г is the boundary of the i-th control unit grid; n is the unit vector in the direction of the external normal corresponding to the boundary Г; S b represents the bottom slope source flux; S f represents the friction source flux; k is the number of the control unit grid edge; l k is the length of the kth side of the i-th control unit grid;

[0031] Step 4.2: Solve the flux term using the HLLC Riemann solver:

[0032] The wave speed is calculated as follows:

[0033]

[0034]

[0035] Among them, the intermediate variable h * and According to the calculation formula,

[0036]

[0037] Step 4.3, the calculation result of F(Q)·n is:

[0038]

[0039] Among them, the control volume unit grid interface intermediate flux F * for,

[0040]

[0041] Control volume unit grid interface flux intermediate wave sides, and for,

[0042]

[0043] In the formula, F L and F R are the fluxes on the left and right sides of the unit grid interface respectively; F * , and are the left wave flux, middle wave flux and right wave flux on both sides of the intermediate wave of the control volume unit grid interface flux respectively; u || is the tangential velocity along the interface, u || =-un y +vn x ,u ||L and u ||R are the left and right tangential velocities of the control volume unit grid interface, in m 3 / s;n x and n y are the indicated directions along x and y, respectively; u and v are the single-width flow rates along x and y, respectively, in m / s; variable Q ⊥ =[h,q x n x +q y n y ];qx and q y are the single width flow along x and y, in m 2 / s;S L , S M , S R They are left wave speed, middle wave speed, and right wave speed, in m / s; u ⊥ is the velocity perpendicular to the grid unit interface, in m / s; u ⊥ =un x +vn y ;h * and is the intermediate variable; g is the acceleration due to gravity, which is 9.81 and the unit is m / s 2 ; C s is the average concentration of the substance along the water depth, kg / m 3 .

[0044] Step 5 is as follows:

[0045] Step 5.1: By introducing the MUSCL variable reconstruction method and Minmod limiter, the numerical discrete format of the convection term of each environmental factor has spatial second-order accuracy and TVD characteristics, as follows:

[0046]

[0047] Where ψ(r) is the Minmod limiter, where r is the pollutant concentration gradient limitation coefficient:

[0048]

[0049] Step 5.2: Treat the diffusion term as a source term and solve it by the second-order central difference method. Thus, the material transport diffusion term is expressed as:

[0050]

[0051] According to the second-order central difference method, the first-order derivative and second-order derivative of the material transport diffusion term are solved, and for any two-dimensional unit grid center (i, j) we get:

[0052]

[0053] Thus, the discrete form of the material transport diffusion term is obtained as follows:

[0054]

[0055] Step 5.3: The bottom slope source term is solved using the bottom slope flux method, and the vector form of the bottom slope flux at a certain side f is obtained as:

[0056]

[0057] In the formula, F sf (Q) is the bottom slope flux at a certain side f; n f is the unit external normal vector of the surface f; n fx and n fy Yes f The component of the unit external normal vector of face f, h L and z bL are the water depth and bottom elevation at the center of the grid cell, respectively. With z ML The water depth and bottom elevation at the edge of the grid. g is the acceleration due to gravity in m / s 2 ;

[0058] Step 5.4: The friction source term is solved using the explicit and implicit friction method. The calculation formula is:

[0059]

[0060] Among them, the calculation formula for single-width flow update to time step n+1 is:

[0061]

[0062] In the formula, S fx and S fy S f Friction source term in x and y directions; C f is the bed surface roughness coefficient, C f =gn 2 / h 1 / 3 , n is the Manning coefficient, unit is s / m 1 / 3 ; h is the water depth, in meters; q x ,q y are the average flow velocities in the x and y directions, in m 2 / s;q x n+1 ,q y n+1 They are the single-width flow rates in the x and y directions at time step n+1, in m 2 / s;q x n ,q y n They are the single-width flow rates in the x and y directions at time step n, in m 2 / s; Δt is the time step, in seconds; Δq x is the difference in single-width flow in the x direction at time step n and time step n+1, that is, Δq x =q x n+1-q x n , unit is m 2 / s.

[0063] Step 6 is as follows:

[0064] Step 6.1: Nitrogen in the nitrogen cycle is considered in both inorganic and organic forms, where organic nitrogen IN is composed of ammonia nitrogen NH3 and nitrate NO3;

[0065] Step 6.2: Comprehensively consider the two forms of phosphorus in the phosphorus cycle: inorganic phosphorus IP and organic phosphorus OP. Inorganic phosphorus IP exists in the water body in the form of phosphate PO4, while organic phosphorus OP is formed by the death of phytoplankton and decomposition of endogenous respiration processes;

[0066] Step 6.3, dissolved oxygen balance includes dissolved oxygen DO, biochemical oxygen demand BOD and chemical oxygen demand COD; dissolved oxygen DO concentration includes oxygen consumption of aquatic plants respiration, nitrification, biodegradation and sediment or silt and oxygen production of photosynthesis, biochemical oxygen demand BOD and chemical oxygen demand COD concentration are related to their own degradation and are degraded with first-order reaction kinetics;

[0067] Step 6.4: Based on the principle of phytoplankton biomass balance, the algae carbon concentration C phlaC represents the phytoplankton biomass and constructs the phytoplankton carbon concentration C phlaC The state equation is

[0068]

[0069] Where, T is water temperature, °C; k cg is the saturated growth rate of phytoplankton at a water temperature of 20°C, d -1 θ cg is the saturated growth temperature coefficient; F(h) is the light limitation function; F(N,P) is the nutrient limitation function in algae photosynthesis.

[0070] Afterwards, the carbon concentration of phytoplankton, C phlaC Relationship between C and chlorophyll α concentration phlaC =k CChla ·C Chla , and obtain the chlorophyll concentration value.

[0071] Step 7 is as follows:

[0072] Step 7.1: Based on the hydrostatic reconstruction method, the bottom slope elevation at the center point M of the dry-wet unit boundary edge f is:

[0073]

[0074] In the formula, are the bottom slope elevations on the left and right sides of the center point M of the dry and wet unit boundary edge f, m;

[0075] Reconstruct the water depth on both sides and make the water level on the right side of the dry-wet interface the same as the elevation of the reconstructed bottom slope, that is:

[0076]

[0077] In the formula, They are the water levels on the left and right sides of the center point M of the dry and wet unit boundary edge f, in m; are the water depths on the left and right sides after reconstruction, in meters; z bM is the bottom slope elevation after reconstruction of the center point M of the dry-wet unit boundary edge f, in meters;

[0078] If expressed in terms of water level, the reconstructed water level is:

[0079]

[0080] In the formula, are the water levels after reconstruction, in m;

[0081] Step 7.2: Based on the above, reconstruct the traffic on both sides of the interface:

[0082]

[0083] In the formula, They are the single-width flows on the left and right sides in the x and y directions after reconstruction, in m 2 / s; are the flow velocities on the left and right sides of the center point M of the dry and wet unit boundary edge f in the x and y directions, respectively, in m / s;

[0084] Step 7.3: To ensure that the reconstructed bottom slope elevation is not higher than the water level on the left, the bottom slope elevation value at the center point M of the dry-wet boundary cell edge f should be the lowest value of the bottom elevation after reconstruction and the water level on the left, that is:

[0085]

[0086] Step 7.4: Further correct u M 、v M , to ensure that the boundary edge of the dry and wet unit changes monotonically, the correction formula is as follows:

[0087]

[0088] In the formula, q xM ,q yM They are the single-width flow rates in the x and y directions after reconstruction of the center point M of the dry and wet unit boundary edge f, in m 3 / s;u M、v M are the corresponding flow velocities in the x and y directions, m / s; h M is the corresponding water depth after reconstruction, in meters; ε is the allowable tolerance, which is given by ε=T 2 It is obtained that T is the area of ​​the largest unit;

[0089] On complex grids, the water depth change value and water depth at the dry-wet boundary interface are used as the discrimination condition to meet the dry-wet boundary solution problem in the second-order or higher-order accuracy format of complex grid terrain. The discrimination formula is as follows:

[0090]

[0091] Where: is the water depth on the left side of the center point M of the dry-wet unit boundary edge f, in meters; is the slope elevation to the left of the center point M of the dry-wet unit boundary edge f, in meters; h i is the water depth of the i-th unit grid, in meters; is the bottom slope elevation on the left side of the i-unit grid, in meters; γ is the ratio of the water depth reconstruction value at the edge of the grid unit to the water depth reconstruction value at the centroid of the unit, γ < 1, the higher the γ value, the more grid units meet this standard, ε wd is the water depth tolerance.

[0092] Step 8 is as follows:

[0093] Step 8.1, using the MUSCL reconstruction method to calculate the spatial variation gradient of each variable in the unit grid; obtaining the restricted variable gradient by restricting the variable gradient, and interpolating the variable value at the midpoint of the unit boundary based on the restricted gradient to achieve MUSCL second-order reconstruction;

[0094] Step 8.2: Obtain the second-order time accuracy by the two-step Runge-Kutta method. The values ​​are:

[0095]

[0096] Intermediate variable value for

[0097]

[0098] Where i is the number of control unit grids; Ω is the volume of the i-th control unit grid, in m 3 ;q n+1 is the single-width flow at time step n+1, in m 2 / s,q n is the single-width flow at time step n, in m 2 / s; Δt is the time step, in seconds; S(q n ) is the source term at time step n; during the entire solution process, GPU parallel computing technology is introduced to improve the computing efficiency;

[0099] Step 8.3: Under the spatial step length Δx, the calculation formula for the value of the time step length Δt is as follows:

[0100]

[0101] In the formula, Δx min It represents the minimum distance from the center of the grid unit to the corresponding interface, in meters; u is the flow velocity, in meters per second; g is the gravitational acceleration, with a value of 9.81 meters per second 2 ; h is the water depth, in meters.

[0102] Step 9 is as follows:

[0103] Step 9.1, specific nutritional status calculation formula:

[0104]

[0105] In the formula, TLI(C chla ) is the chlorophyll α nutritional status index of phytoplankton; TLI (C TP ) represents the total phosphorus nutritional status index; TLI (C TN ) is the total nitrogen nutritional status index; TLI (C COD ) is the chemical oxygen demand nutrient status index; C chla is the concentration of chlorophyll α, in mg / L; C TP is the total phosphorus concentration, in mg / L; C TN is the total nitrogen concentration, in mg / L; C COD is the chemical oxygen demand concentration, in mg / L;

[0106] Step 9.2: Using the above calculation formulas for each nutritional status index and considering the weight of each nutritional status index, the formula for the comprehensive nutritional status index is:

[0107]

[0108] Where, TLI(∑) is the comprehensive nutrient status index determined by the above environmental factors, reflecting the eutrophication of rivers and lakes; TLI(j) is the nutrient status index of the jth parameter; W j is the corresponding weight value of the j-th parameter nutritional status index;

[0109] Taking the chlorophyll α concentration of phytoplankton as the benchmark, the relevant weight of the jth parameter is obtained, and the normalized formula is:

[0110]

[0111] In the formula, is the correlation coefficient between the jth parameter and the chlorophyll α concentration; m is the number of parameters, which is used to evaluate the eutrophication of river and lake water bodies.

[0112] The beneficial effect of the present invention is that it provides an efficient water eutrophication assessment method that integrates inorganic and organic substances. In view of the disadvantage that the existing water eutrophication assessment model involves many ecosystem state variables and the mutual conversion relationship between them is relatively complex, resulting in low model simulation efficiency, the proposed efficient water eutrophication assessment method that comprehensively considers inorganic and organic substances is aimed at comprehensively considering the organic and inorganic forms of nitrogen and phosphorus substances in water bodies, mastering the material, chemical and biological processes inside water bodies to explain the occurrence mechanism of river and lake eutrophication, so as to evaluate the effectiveness of the implementation of engineering measures in changing the eutrophic state of water bodies and the corresponding water quality improvement. BRIEF DESCRIPTION OF THE DRAWINGS

[0113] Figure 1 It is a schematic diagram of the study area and its high-resolution DEM topographic map;

[0114] Figure 2 It is a schematic diagram of an efficient water eutrophication model that comprehensively considers inorganic and organic substances;

[0115] Figure 3 This is a comparison chart of the simulated and measured values ​​of various environmental factor variables at the Chanba River monitoring points;

[0116] Figure 4 It is a quantitative evaluation result diagram of the performance of the high-efficiency water eutrophication assessment model that comprehensively considers inorganic and organic substances;

[0117] Figure 5(a) shows the flow rate of 7.15m 3 / s condition, the water depth and flow velocity distribution diagram of Chanba River;

[0118] Figure 5(b) shows the flow rate of 7.15m 3 / s conditions. DETAILED DESCRIPTION

[0119] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0120] The present invention provides a highly efficient water eutrophication assessment method that integrates inorganic and organic substances, and is specifically implemented according to the following steps:

[0121] Step 1: Using ArcGIS software to obtain DEM terrain data as the input condition of the two-dimensional hydrodynamic-material transport coupling model, so as to characterize the real terrain conditions of rivers, lakes and seas;

[0122] Step 2: Using the rectangular grid discretized DEM terrain data, assign values ​​to the relevant parameters in the model according to the corresponding known boundary conditions;

[0123] Step 3: Based on the material convection diffusion theory and the principle of mutual biochemical reaction between materials, construct the material transport diffusion equation;

[0124] The material transport and diffusion equation in step 3 is expressed as follows:

[0125]

[0126] In the formula, C s is the average substance concentration along the water depth direction, in kg / m 3 ;D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s; S in source and sink items sc h is the biochemical reaction term, and its unit is kg / m 3 ; C0 is the initial average substance concentration along the water depth direction, in kg / m 3 .

[0127] Step 4: Couple the material transport equation with the high-precision hydrodynamic model to construct a two-dimensional hydrodynamic-material transport coupling control equation group, which is expressed in vector form; the model flux is discretized using the central finite volume method in the Godunov format, which includes the HLLC approximate solution method for integrated material convection and the second-order central difference method for its diffusion term to overcome the numerical oscillation and dissipation caused by solving material convection, so as to achieve synchronization of the coupled solution of hydrodynamics and material transport.

[0128] The vector form of the two-dimensional hydrodynamic-mass transport coupling equations in step 4 is as follows:

[0129]

[0130] in,

[0131]

[0132] Where t is time in seconds; Q is a variable vector, where h is water depth in meters, and q is x and q y They are the single width flow in the x and y directions, in m 2 / s; F, G are flux vectors in the x and y directions respectively; u, v are flow velocities in the x and y directions respectively, in units of m / s; S is the source term vector, including bottom slope source term, friction source term, material diffusion term and material source term; zb is the riverbed bottom elevation, in m; C f is the friction coefficient of the bed surface, calculated as C f =gn 2 / h 1 / 3 , where n is the Manning coefficient; g is the acceleration due to gravity, in m / s 2 ; C s is the average substance concentration along the water depth direction, in kg / m 3 ;D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s; S in source and sink items sc h is the biochemical reaction term, and its unit is kg / m 3 ; C0 is the initial average substance concentration along the water depth direction, in kg / m 3 .

[0133] Step 4 is as follows:

[0134] Step 4.1: In the control unit grid i, the central finite volume method is used to integrate the hydrodynamic-material transport diffusion coupling control equations. Based on the Gaussian divergence theorem, the area integral of the flux term is converted into a line integral in the form of:

[0135]

[0136] Where Ω is the volume of control volume i; Г is the boundary of control volume i; n is the unit vector in the direction of the external normal corresponding to the control volume boundary Г; F(Q)·n is the external normal convection term flux vector on the control volume boundary Г; S is the source term vector, including the bottom slope source term, friction source term, diffusion term and biochemical reaction term,

[0137] The flux F(Q)·n line integral is:

[0138]

[0139] Where Ω is the volume of the i-th control unit grid, in m 3 ; Г is the boundary of the i-th control unit grid; n is the unit vector in the direction of the external normal corresponding to the boundary Г; S b represents the bottom slope source flux; S f represents the friction source flux; k is the number of the control unit grid edge; l k is the length of the kth side of the i-th control unit grid;

[0140] Step 4.2: Solve the flux term using the HLLC Riemann solver:

[0141] The wave speed is calculated as follows:

[0142]

[0143]

[0144] Among them, the intermediate variable h * and According to the calculation formula,

[0145]

[0146] Step 4.3, the calculation result of F(Q)·n is:

[0147]

[0148] Among them, the control volume unit grid interface intermediate flux F * for,

[0149]

[0150] Control volume unit grid interface flux intermediate wave sides, and for,

[0151]

[0152] In the formula, F L and F R are the fluxes on the left and right sides of the unit grid interface respectively; F * , and are the left wave flux, middle wave flux and right wave flux on both sides of the intermediate wave of the control volume unit grid interface flux respectively; u || is the tangential velocity along the interface, u || =-un y +vn x ,u ||L and u ||R are the left and right tangential velocities of the control volume unit grid interface, in m 3 / s;n x and n y are the indicated directions along x and y, respectively; u and v are the single-width flow rates along x and y, respectively, in m / s; variable Q ⊥ =[h,q x n x +q y n y ];q x and q y are the single width flow along x and y, in m 2 / s;S L , SM , S R They are left wave speed, middle wave speed, and right wave speed, in m / s; u ⊥ is the velocity perpendicular to the grid unit interface, in m / s; u ⊥ =un x +vn y ;h * and is the intermediate variable; g is the acceleration due to gravity, which is 9.81 and the unit is m / s 2 ; C s is the average concentration of the substance along the water depth, kg / m 3 .

[0153] Step 5: In order to make the numerical discrete format of water quality convection terms have spatial second-order accuracy and TVD characteristics, the MUSCL variable reconstruction method and Minmod limiter are introduced in the calculation of water quality Riemann variables to limit the concentration gradient of each environmental factor in the model. In addition, other terms in the equation group are further solved;

[0154] Step 5 is as follows:

[0155] Step 5.1: By introducing the MUSCL variable reconstruction method and Minmod limiter, the numerical discrete format of the convection term of each environmental factor has spatial second-order accuracy and TVD characteristics, as follows:

[0156]

[0157] Where ψ(r) is the Minmod limiter, where r is the pollutant concentration gradient limitation coefficient:

[0158]

[0159] Step 5.2: Treat the diffusion term as a source term and solve it by the second-order central difference method. Thus, the material transport diffusion term is expressed as:

[0160]

[0161] According to the second-order central difference method, the first-order derivative and second-order derivative of the material transport diffusion term are solved, and for any two-dimensional unit grid center (i, j) we get:

[0162]

[0163] Thus, the discrete form of the material transport diffusion term is obtained as follows:

[0164]

[0165] Step 5.3: The bottom slope source term is solved using the bottom slope flux method, and the vector form of the bottom slope flux at a certain side f is obtained as:

[0166]

[0167] In the formula, F sf (Q) is the bottom slope flux at a certain side f; n f is the unit external normal vector of the surface f; n fx and n fy Yes f The component of the unit external normal vector of face f, h L and z bL are the water depth and bottom elevation at the center of the grid cell, respectively. With z ML The water depth and bottom elevation at the edge of the grid. g is the acceleration due to gravity in m / s 2 ;

[0168] Step 5.4: The friction source term is solved using the explicit and implicit friction method. The calculation formula is:

[0169]

[0170] Among them, the calculation formula for single-width flow update to time step n+1 is:

[0171]

[0172]

[0173] In the formula, S fx and S fy S f Friction source term in x and y directions; C f is the bed surface roughness coefficient, C f =gn 2 / h 1 / 3 , n is the Manning coefficient, unit is s / m 1 / 3 ; h is the water depth, in meters; q x ,q y are the average flow velocities in the x and y directions, in m 2 / s;q x n+1 ,q y n+1 They are the single-width flow rates in the x and y directions at time step n+1, in m 2 / s;q x n ,q y n They are the single-width flow rates in the x and y directions at time step n, in m2 / s; Δt is the time step, in seconds; Δq x is the difference in single-width flow in the x direction at time step n and time step n+1, that is, Δq x =q x n+1 -q x n , unit is m 2 / s; In order to avoid the denominator being zero, θ is taken as 1.0e -12 .

[0174] Step 6: Substitute the biochemical reaction term S in the source term of the hydrodynamic-material transport coupling equation sc h, and the dissection operator method is used to convert it into the balance equation of each material state variable and solve it one by one. In this balance equation, the inorganic and organic forms of nutrients such as nitrogen and phosphorus in the water body are comprehensively considered. That is, inorganic nitrogen (ammonia nitrogen NH3, nitrate nitrogen NO3) and organic nitrogen (ON) are considered in the nitrogen cycle; inorganic phosphorus (phosphate PO4) and organic phosphorus (OP) are considered in the phosphorus cycle. At the same time, material state variables such as biochemical oxygen demand BOD, chemical oxygen demand COD, dissolved oxygen DO, and algae carbon concentration are considered.

[0175] Step 6 is as follows:

[0176] Step 6.1: Nitrogen in the nitrogen cycle is considered in both inorganic and organic forms, where organic nitrogen IN is composed of ammonia nitrogen NH3 and nitrate NO3;

[0177] Step 6.2: Inorganic phosphorus IP and organic phosphorus OP are considered in the phosphorus cycle. Inorganic phosphorus IP only exists in the form of phosphate PO4 in the water body, while organic phosphorus OP is formed by the death of phytoplankton and decomposition of endogenous respiration process.

[0178] Step 6.3, dissolved oxygen balance includes dissolved oxygen DO, biochemical oxygen demand BOD and chemical oxygen demand COD; dissolved oxygen DO concentration includes oxygen consumption of aquatic plants respiration, nitrification, biodegradation and sediment or silt and oxygen production of photosynthesis, biochemical oxygen demand BOD and chemical oxygen demand COD concentration are related to their own degradation and are degraded with first-order reaction kinetics;

[0179] Step 6.4: Based on the principle of phytoplankton biomass balance, the algae carbon concentration C phlaC represents the phytoplankton biomass and constructs the phytoplankton carbon concentration C phlaC The state equation is

[0180]

[0181] Where, T is water temperature, °C; k cgis the saturated growth rate of phytoplankton at a water temperature of 20°C, d -1 θ cg is the saturated growth temperature coefficient; F(h) is the light limitation function; F(N,P) is the nutrient limitation function in algae photosynthesis.

[0182] Afterwards, the carbon concentration of phytoplankton, C phlaC Relationship between C and chlorophyll α concentration phlaC =k CChla ·C Chla , and obtain the chlorophyll concentration value.

[0183] Step 7: The bottom slope elevation of the unit grid boundary surface is defined as the maximum value at the boundary between two adjacent units through the improved hydrostatic reconstruction method, and the corresponding water depth and water level are reconstructed. To maintain the negative water depth value, the negative water depth is defined as zero to maintain mass conservation and stability at the dry-wet boundary.

[0184] Step 7 is as follows:

[0185] Step 7.1: Based on the hydrostatic reconstruction method, the bottom slope elevation at the center point M of the dry-wet unit boundary edge f is: In the formula, are the bottom slope elevations on the left and right sides of the center point M at the boundary edge f of the dry and wet unit, respectively;

[0186] Reconstruct the water depth on both sides and make the water level on the right side of the dry-wet interface the same as the elevation of the reconstructed bottom slope, that is:

[0187]

[0188] In the formula, They are the water levels on the left and right sides of the center point M of the dry and wet unit boundary edge f, in m; are the water depths on the left and right sides after reconstruction, in meters; z bM is the bottom slope elevation after reconstruction of the center point M of the dry-wet unit boundary edge f, in meters;

[0189] If expressed in terms of water level, the reconstructed water level is:

[0190]

[0191] In the formula, are the water levels after reconstruction, in m;

[0192] Step 7.2: Based on the above, reconstruct the traffic on both sides of the interface:

[0193]

[0194] In the formula, They are the single-width flows on the left and right sides in the x and y directions after reconstruction, in m 2 / s; are the flow velocities on the left and right sides of the center point M of the dry and wet unit boundary edge f in the x and y directions, respectively, in m / s;

[0195] Step 7.3: To ensure that the reconstructed bottom slope elevation is not higher than the water level on the left, the bottom slope elevation value at the center point M of the dry-wet boundary cell edge f should be the lowest value of the bottom elevation after reconstruction and the water level on the left, that is:

[0196]

[0197] Step 7.4: Further correct u M 、v M , to ensure that the boundary edge of the dry and wet unit changes monotonically, the correction formula is as follows:

[0198]

[0199] In the formula, q xM ,q yM They are the single-width flow rates in the x and y directions after reconstruction of the center point M of the dry and wet unit boundary edge f, in m 3 / s;u M 、v M are the corresponding flow velocities in the x and y directions, m / s; h M is the corresponding water depth after reconstruction, in meters; ε is the allowable tolerance, which is given by ε=T 2 It is obtained that T is the area of ​​the largest unit;

[0200] On complex grids, the water depth change value and water depth at the dry-wet boundary interface are used as the discrimination condition to meet the dry-wet boundary solution problem in the second-order or higher-order accuracy format of complex grid terrain. The discrimination formula is as follows:

[0201]

[0202] Where: is the water depth on the left side of the center point M of the dry-wet unit boundary edge f, in meters; is the slope elevation to the left of the center point M of the dry-wet unit boundary edge f, in meters; h i is the water depth of the i-th unit grid, in meters; is the bottom slope elevation on the left side of the i-unit grid, in meters; γ is the ratio of the water depth reconstruction value at the edge of the grid unit to the water depth reconstruction value at the centroid of the unit, γ < 1, the higher the γ value, the more grid units meet this standard, and it is taken as 0.25 in this paper; ε wd is the water depth tolerance, usually ε wd =1×10 -6 .

[0203] Step 8: Based on the above, the second-order precision of MUSCL and the second-order Runge-Kutta method are further used to achieve the second-order spatiotemporal precision of the model, and the CFL condition is used to restrict the temporal stability. At the same time, GPU parallel computing is introduced in the entire model simulation process to improve the calculation speed of the model.

[0204] Step 8 is as follows:

[0205] Step 8.1, using the MUSCL reconstruction method to calculate the spatial variation gradient of each variable in the unit grid; obtaining the restricted variable gradient by restricting the variable gradient, and interpolating the variable value at the midpoint of the unit boundary based on the restricted gradient to achieve MUSCL second-order reconstruction;

[0206] Step 8.2: Obtain the second-order time accuracy by the two-step Runge-Kutta method. The values ​​are:

[0207]

[0208] Intermediate variable value for

[0209]

[0210] Where i is the number of control unit grids; Ω is the volume of the i-th control unit grid, in m 3 ;q n+1 is the single-width flow at time step n+1, in m 2 / s,q n is the single-width flow at time step n, in m 2 / s; Δt is the time step, in seconds; S(q n ) is the source term at time step n; during the entire solution process, GPU parallel computing technology is introduced to improve the computing efficiency;

[0211] Step 8.3: Under the spatial step length Δx, the calculation formula for the value of the time step length Δt is as follows:

[0212]

[0213] In the formula, Δx min It represents the minimum distance from the center of the grid unit to the corresponding interface, in meters; u is the flow velocity, in meters per second; g is the gravitational acceleration, with a value of 9.81 meters per second 2 ; h is the water depth, in meters.

[0214] Step 9: Based on the calculated concentrations of various environmental factors, taking the chlorophyll α (CHLA) concentration of phytoplankton as a benchmark, and considering total nitrogen, total phosphorus and chemical oxygen demand (COD), a comprehensive nutrient status index method is constructed to evaluate the eutrophication of river and lake water bodies;

[0215] Step 9 is as follows:

[0216] Step 9.1, specific nutritional status calculation formula:

[0217]

[0218] In the formula, TLI(C chla ) is the chlorophyll α nutritional status index of phytoplankton; TLI (C TP ) represents the total phosphorus nutritional status index; TLI (C TN ) is the total nitrogen nutritional status index; TLI (C COD ) is the chemical oxygen demand nutrient status index; C chla is the concentration of chlorophyll α, in mg / L; C TP is the total phosphorus concentration, in mg / L; C TN is the total nitrogen concentration, in mg / L; C COD is the chemical oxygen demand concentration, in mg / L;

[0219] Step 9.2: Using the above calculation formulas for each nutritional status index and considering the weight of each nutritional status index, the formula for the comprehensive nutritional status index is:

[0220]

[0221] Where, TLI(Σ) is the comprehensive nutrient status index determined by the above environmental factors, reflecting the eutrophication of rivers and lakes; TLI(j) is the nutrient status index of the jth parameter; W j is the corresponding weight value of the j-th parameter nutritional status index;

[0222] Taking the chlorophyll α concentration of phytoplankton as the benchmark, the relevant weight of the jth parameter is obtained, and the normalized formula is:

[0223] In the formula, is the correlation coefficient between the jth parameter and the chlorophyll α concentration; m is the number of parameters, where The selection is based on the results of the lake survey in my country, in order to evaluate the eutrophication of river and lake water bodies.

[0224] Step 10: Based on the concentration change process of each environmental factor obtained by the above model simulation, the MATLAB program is used to output the distribution diagram of the transport process of each environmental factor to intuitively reflect the transport law of each environmental factor in the water body, so as to explain in detail the occurrence mechanism of river and lake eutrophication. At the same time, the comprehensive index method is used to evaluate the eutrophication status of the water body.

[0225] Example 1

[0226] The present invention provides a highly efficient water eutrophication assessment method that integrates inorganic and organic substances, and is specifically implemented according to the following steps:

[0227] Step 1: Using ArcGIS software to obtain DEM terrain data as the input condition of the two-dimensional hydrodynamic-material transport coupling model, so as to characterize the real terrain conditions of rivers, lakes and seas;

[0228] Step 2: Using the rectangular grid discretized DEM terrain data, assign values ​​to the relevant parameters in the model according to the corresponding known boundary conditions;

[0229] Step 3: Construct the material transport and diffusion equation;

[0230] Step 4: Couple the material transport equation with the high-precision hydrodynamic model to construct a two-dimensional hydrodynamic-material transport coupling control equation set.

[0231] Step 5: Introduce the MUSCL variable reconstruction method and Minmod limiter in the calculation of water quality Riemann variables to limit the concentration gradient of each environmental factor in the model;

[0232] Step 6: Substitute the biochemical reaction term S in the source term of the hydrodynamic-material transport coupling equation sc h, and use the split operator method to transform it into the transformation equilibrium equation of each material state variable and solve it one by one.

[0233] Step 7: The bottom slope elevation of the unit grid boundary surface is defined as the maximum value at the boundary between two adjacent units through the improved hydrostatic reconstruction method, and the corresponding water depth and water level are reconstructed.

[0234] Step 8: Use the MUSCL second-order accuracy and the second-order Runge-Kutta method to achieve the second-order spatiotemporal accuracy of the model.

[0235] Step 9: Based on the calculated concentrations of various environmental factors and taking the chlorophyll α (CHLA) concentration of phytoplankton as a benchmark, a comprehensive nutrient status index method is constructed to evaluate the eutrophication of river and lake water bodies;

[0236] Step 10: Output the distribution diagram of the transport process of each environmental factor to intuitively reflect the transport law of each environmental factor in the water body, and evaluate the eutrophication status of the water body through the comprehensive index method.

[0237] Example 2

[0238] The present invention provides a highly efficient water eutrophication assessment method that integrates inorganic and organic substances, and is specifically implemented according to the following steps:

[0239] Step 1: Using ArcGIS software to obtain DEM terrain data as the input condition of the two-dimensional hydrodynamic-material transport coupling model, so as to characterize the real terrain conditions of rivers, lakes and seas;

[0240] Step 2: Using the rectangular grid discretized DEM terrain data, assign values ​​to the relevant parameters in the model according to the corresponding known boundary conditions;

[0241] Step 3: Construct the material transport and diffusion equation;

[0242] Step 4: Couple the material transport equation with the high-precision hydrodynamic model to construct a two-dimensional hydrodynamic-material transport coupling control equation set.

[0243] Step 5: Introduce the MUSCL variable reconstruction method and Minmod limiter in the calculation of water quality Riemann variables to limit the concentration gradient of each environmental factor in the model;

[0244] Step 6: Substitute the biochemical reaction term S in the source term of the hydrodynamic-material transport coupling equation sc h, and use the split operator method to transform it into the transformation equilibrium equation of each material state variable and solve it one by one.

[0245] Step 7: The bottom slope elevation of the unit grid boundary surface is defined as the maximum value at the boundary between two adjacent units through the improved hydrostatic reconstruction method, and the corresponding water depth and water level are reconstructed.

[0246] Step 8: Based on the above, the second-order accuracy of MUSCL and the second-order Runge-Kutta method are further used to achieve the second-order spatiotemporal accuracy of the model;

[0247] Step 8 is as follows:

[0248] Step 8.1, using the MUSCL reconstruction method to calculate the spatial variation gradient of each variable in the unit grid; obtaining the restricted variable gradient by restricting the variable gradient, and interpolating the variable value at the midpoint of the unit boundary based on the restricted gradient to achieve MUSCL second-order reconstruction;

[0249] Step 8.2: Obtain the second-order time accuracy by the two-step Runge-Kutta method. The values ​​are:

[0250] Intermediate variable value for:

[0251]

[0252] Where i is the number of control unit grids; Ω is the volume of the i-th control unit grid, in m 3 ;q n+1 is the single-width flow at time step n+1, in m 2 / s,q n is the single-width flow at time step n, in m 2 / s; Δt is the time step, in seconds; S(q n ) is the source term at time step n; during the entire solution process, GPU parallel computing technology is introduced to improve the computing efficiency;

[0253] Step 8.3: Under the spatial step length Δx, the calculation formula for the value of the time step length Δt is as follows:

[0254]

[0255] In the formula, Δx min It represents the minimum distance from the center of the grid unit to the corresponding interface, in meters; u is the flow velocity, in meters per second; g is the gravitational acceleration, with a value of 9.81 meters per second 2 ; h is the water depth, in meters.

[0256] Step 9: Based on the calculated concentrations of various environmental factors and taking the chlorophyll α (CHLA) concentration of phytoplankton as a benchmark, a comprehensive nutrient status index method is constructed to evaluate the eutrophication of river and lake water bodies;

[0257] Step 10: Output the distribution diagram of the transport process of each environmental factor to intuitively reflect the transport law of each environmental factor in the water body, and evaluate the eutrophication status of the water body through the comprehensive index method.

[0258] Example 3

[0259] This method takes the Xi'an Chanba River Basin as an example for model application. The Xi'an Chanba River Basin is located in the southeast of Xi'an (109°00'E~109°47'E, 33°50"N~34°27"N), about 78km long from north to south, about 50km wide from east to west, and has a total area of ​​2581km 2 The comprehensive roughness of the river channel is 0.035s / m 1 / 3 Among them, the Ba River is a first-level tributary of the Wei River. On average, its runoff is 743 million m 3 , and the average gradient is 6.2‰. The river originates from Jiudaogou in Bayuan Township, Lantian County, with a total length of 104.1km. It flows to Sanlang Village in Baqiao District and flows into the Weihe River. The basin area of ​​the tributary of the river, Chanhe River, is 760km 2, the average runoff over many years is 188 million m 3 The total length is 64 km and the average gradient is 8.9‰. It joins the river at Tanjiabao in Weiyang District. The location of Chanba Ecological Zone and the topographic map of the studied river section can be found at Figure 1 The terrain of Yanming Lake No. 5 was discretized using rectangular grids with a grid size of 10.0 m, totaling 683,928 grids.

[0260] Combination Figure 2 In the simulation, the measured data of water environmental factors (DO, NH3-N, COD, TN, TP and chlorophyll α) at ​​the monitoring point of Chanba River were selected to calibrate the relevant parameters. The comparison chart of verification results is shown in Figure 3 As can be seen from the figure, the concentration change trend of each water environmental factor is basically consistent with the measured data. The measured values ​​are distributed near the simulation curve, and the two are in good agreement. At the same time, in order to illustrate the simulation effect of the model, based on the measured and simulated values ​​of the state variables at the above monitoring points, combined with Figure 4 Figure 5(a) and Figure 5(b) use three statistical variables, MAE, RMSE and NSE, to quantitatively evaluate the simulation results. According to the evaluation results, the error statistical variables MAE, RMSE and NSE of each water environment factor are 0.002-0.943, 0.005-3.001 and 0.740-0.987 respectively. Combined with the above comparison charts, it can be seen that the simulation accuracy of the model is good and can truly reflect the water environment changes in the Chanba River. Based on the calibration of these relevant parameters, a flow of 7.5m 3 / s conditions, the distribution of hydrodynamics and environmental factors of the Chanba River, in order to reflect the distribution law of environmental factor variables in the Chanba River.

Claims

1. An efficient water eutrophication assessment method integrating inorganic and organic substances, characterized in that: Follow the steps below to implement it: Step 1: Obtain DEM terrain data as input conditions for the two-dimensional hydrodynamic-material transport coupling model; Step 2: Use the rectangular grid discretized DEM terrain data to assign values ​​to the relevant parameters in the model; Step 3: Construct the material transport and diffusion equation; Step 4: construct a two-dimensional hydrodynamic-material transport coupling control equation set; Step 5: Limit the concentration gradient of each environmental factor in the model; Step 6: Substitute the biochemical reaction term S in the source term of the hydrodynamic-material transport coupling equation sc h is converted into the state variable transformation equilibrium equation of each substance and solved one by one; Step 7: define the bottom slope elevation of the unit grid boundary surface as the maximum value at the boundaries of two adjacent units for reconstruction, and reconstruct the corresponding water depth and water level; Step 8: Achieve the spatiotemporal second-order accuracy of the model; Step 9: Evaluate the eutrophication of rivers and lakes; Step 10: Evaluate the eutrophication status of water bodies through the comprehensive index method.

2. The efficient water eutrophication assessment method of integrated inorganic and organic substances according to claim 1 is characterized in that: The material transport diffusion equation in step 3 is expressed as follows: In the formula, C s is the average substance concentration along the water depth direction, in kg / m 3 ; D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s; S in source and sink items sc h is the biochemical reaction term, and its unit is kg / m 3 ; C0 is the initial average substance concentration along the water depth direction, in kg / m 3 .

3. The efficient water eutrophication assessment method of integrated inorganic and organic substances according to claim 2 is characterized in that: The vector form of the two-dimensional hydrodynamic-material transport coupling equations in step 4 is expressed as follows: in, Where t is time in seconds; Q is a variable vector, where h is water depth in meters, and q is x and q y They are the single width flow in the x and y directions, in m 2 / s; F, G are flux vectors in the x and y directions respectively; u, v are flow velocities in the x and y directions respectively, in units of m / s; S is the source term vector, including bottom slope source term, friction source term, material diffusion term and material source term; z b is the riverbed bottom elevation, in m; C f is the friction coefficient of the bed surface, calculated as C f =gn 2 / h 1 / 3 , where n is the Manning coefficient; g is the acceleration due to gravity, in m / s 2 ; C s is the average substance concentration along the water depth direction, in kg / m 3 ;D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s; S in source and sink items sc h is the biochemical reaction term, and its unit is kg / m 3 ; C0 is the initial average substance concentration along the water depth direction, in kg / m 3 .

4. The efficient water eutrophication assessment method of integrated inorganic and organic substances according to claim 3 is characterized in that: The step 4 is specifically as follows: Step 4.1: In the control unit grid i, the central finite volume method is used to integrate the hydrodynamic-material transport diffusion coupling control equations. Based on the Gaussian divergence theorem, the area integral of the flux term is converted into a line integral in the form of: Where Ω is the volume of control volume i; Г is the boundary of control volume i; n is the unit vector in the direction of the external normal corresponding to the control volume boundary Г; F(Q)·n is the external normal convection term flux vector on the control volume boundary Г; S is the source term vector, including the bottom slope source term, friction source term, diffusion term and biochemical reaction term, The flux F(Q)·n line integral is: Where Ω is the volume of the i-th control unit grid, in m 3 ; Г is the boundary of the i-th control unit grid; n is the unit vector in the direction of the external normal corresponding to the boundary Г; S b represents the bottom slope source flux; S f represents the friction source flux; k is the number of the control unit grid edge; l k is the length of the kth side of the i-th control unit grid; Step 4.2: Solve the flux term using the HLLC Riemann solver: The wave speed is calculated as follows: Among them, the intermediate variable h * and According to the calculation formula, Step 4.3, the calculation result of F(Q)·n is: Among them, the control volume unit grid interface intermediate flux F * for, Control volume unit grid interface flux intermediate wave sides, and for, In the formula, F L and F R are the fluxes on the left and right sides of the unit grid interface respectively; F * , and are the left wave flux, middle wave flux and right wave flux on both sides of the intermediate wave of the control volume unit grid interface flux respectively; u || is the tangential velocity along the interface, u || =-un y +vn x ,u ||L and u ||R are the left and right tangential velocities of the control volume unit grid interface, in m 3 / s;n x and n y are the indicated directions along x and y, respectively; u and v are the single-width flow rates along x and y, respectively, in m / s; variable Q ⊥ =[h,q x n x +q y n y ];q x and q y are the single width flow along x and y, in m 2 / s;S L , S M , S R They are left wave speed, middle wave speed, and right wave speed, in m / s; u ⊥ is the velocity perpendicular to the grid unit interface, in m / s; u ⊥ =un x +vn y ;h * and is the intermediate variable; g is the acceleration due to gravity, which is 9.81 and the unit is m / s 2 ; C s is the average concentration of the substance along the water depth, kg / m 3 .

5. The efficient water eutrophication assessment method of integrated inorganic and organic substances according to claim 4 is characterized in that: The step 5 is specifically as follows: Step 5.1: By introducing the MUSCL variable reconstruction method and Minmod limiter, the numerical discrete format of the convection term of each environmental factor has spatial second-order accuracy and TVD characteristics, as follows: Where ψ(r) is the Minmod limiter, where r is the pollutant concentration gradient limitation coefficient: Step 5.2: Treat the diffusion term as a source term and solve it by the second-order central difference method. Thus, the material transport diffusion term is expressed as: According to the second-order central difference method, the first-order derivative and second-order derivative of the material transport diffusion term are solved, and for any two-dimensional unit grid center (i, j) we get: Thus, the discrete form of the material transport diffusion term is obtained as follows: Step 5.3: The bottom slope source term is solved using the bottom slope flux method, and the vector form of the bottom slope flux at a certain side f is obtained as: In the formula, F sf (Q) is the bottom slope flux at a certain side f; n f is the unit external normal vector of the surface f; n fx and n fy Yes f The component of the unit external normal vector of face f, h L and z bL are the water depth and bottom elevation at the center of the grid cell, respectively. With z ML The water depth and bottom elevation at the edge of the grid. g is the acceleration due to gravity in m / s 2 ; Step 5.4: The friction source term is solved using the explicit and implicit friction method. The calculation formula is: Among them, the calculation formula for single-width flow update to time step n+1 is: In the formula, S fx and S fy S f Friction source term in x and y directions; C f is the bed surface roughness coefficient, C f =gn 2 / h 1 / 3 , n is the Manning coefficient, unit is s / m 1 / 3 ; h is the water depth, in meters; q x ,q y are the average flow velocities in the x and y directions, in m 2 / s;q x n +1 ,q y n+1 They are the single-width flow rates in the x and y directions at time step n+1, in m 2 / s;q x n ,q y n They are the single-width flow rates in the x and y directions at time step n, in m 2 / s; Δt is the time step, in seconds; Δq x is the difference in single-width flow in the x direction at time step n and time step n+1, that is, Δq x =q x n+1 -q x n , unit is m 2 / s.

6. The efficient water eutrophication assessment method of integrated inorganic and organic substances according to claim 5 is characterized in that: The step 6 is specifically as follows: Step 6.1: Nitrogen in the nitrogen cycle is considered in both inorganic and organic forms, where organic nitrogen IN is composed of ammonia nitrogen NH3 and nitrate NO3; Step 6.2: Comprehensively consider the two forms of phosphorus in the phosphorus cycle: inorganic phosphorus IP and organic phosphorus OP. Inorganic phosphorus IP exists in the water body in the form of phosphate PO4, while organic phosphorus OP is formed by the death of phytoplankton and decomposition of endogenous respiration processes; Step 6.3, dissolved oxygen balance includes dissolved oxygen DO, biochemical oxygen demand BOD and chemical oxygen demand COD; dissolved oxygen DO concentration includes oxygen consumption of aquatic plants respiration, nitrification, biodegradation and sediment or silt and oxygen production of photosynthesis, biochemical oxygen demand BOD and chemical oxygen demand COD concentration are related to their own degradation and are degraded with first-order reaction kinetics; Step 6.4: Based on the principle of phytoplankton biomass balance, the algae carbon concentration C phlaC represents the phytoplankton biomass and constructs the phytoplankton carbon concentration C phlaC The state equation is Where, T is water temperature, °C; k cg is the saturated growth rate of phytoplankton at a water temperature of 20°C, d -1 θ cg is the saturated growth temperature coefficient; F(h) is the light limitation function; F(N,P) is the nutrient limitation function in algae photosynthesis; Afterwards, the carbon concentration of phytoplankton, C phlaC Relationship between C and chlorophyll α concentration phlaC =k CChla ·C Chla , and obtain the chlorophyll concentration value.

7. The efficient water eutrophication assessment method of integrated inorganic and organic substances according to claim 6 is characterized in that: The step 7 is specifically as follows: Step 7.1: Based on the hydrostatic reconstruction method, the bottom slope elevation at the center point M of the dry-wet unit boundary edge f is: In the formula, are the bottom slope elevations on the left and right sides of the center point M of the dry and wet unit boundary edge f, m; Reconstruct the water depth on both sides and make the water level on the right side of the dry-wet interface the same as the elevation of the reconstructed bottom slope, that is: In the formula, They are the water levels on the left and right sides of the center point M of the dry and wet unit boundary edge f, in m; are the water depths on the left and right sides after reconstruction, in meters; z bM is the bottom slope elevation after reconstruction of the center point M of the dry-wet unit boundary edge f, in meters; If expressed in terms of water level, the reconstructed water level is: In the formula, are the water levels after reconstruction, in m; Step 7.2: Based on the above, reconstruct the traffic on both sides of the interface: In the formula, They are the single-width flows on the left and right sides in the x and y directions after reconstruction, in m 2 / s; are the flow velocities on the left and right sides of the center point M of the dry and wet unit boundary edge f in the x and y directions, respectively, in m / s; Step 7.3: To ensure that the reconstructed bottom slope elevation is not higher than the water level on the left, the bottom slope elevation value at the center point M of the dry-wet boundary cell edge f should be the lowest value of the bottom elevation after reconstruction and the water level on the left, that is: Step 7.4: Further correct u M 、v M , to ensure that the boundary edge of the dry and wet unit changes monotonically, the correction formula is as follows: In the formula, q xM ,q yM They are the single-width flow rates in the x and y directions after reconstruction of the center point M of the dry and wet unit boundary edge f, in m 3 / s;u M 、v M are the corresponding flow velocities in the x and y directions, m / s; h M is the corresponding water depth after reconstruction, in meters; ε is the allowable tolerance, which is given by ε=T 2 It is obtained that T is the area of ​​the largest unit; On complex grids, the water depth change value and water depth at the dry-wet boundary interface are used as the discrimination condition to meet the dry-wet boundary solution problem in the second-order or higher-order accuracy format of complex grid terrain. The discrimination formula is as follows: Where: is the water depth on the left side of the center point M of the dry-wet unit boundary edge f, in meters; is the slope elevation to the left of the center point M of the dry-wet unit boundary edge f, in meters; h i is the water depth of the i-th unit grid, in meters; is the bottom slope elevation on the left side of the i-unit grid, in meters; γ is the ratio of the water depth reconstruction value at the edge of the grid unit to the water depth reconstruction value at the centroid of the unit, γ < 1, the higher the γ value, the more grid units meet this standard, ε wd is the water depth tolerance.

8. The efficient water eutrophication assessment method based on comprehensive inorganic and organic substances according to claim 7 is characterized in that: The step 8 is specifically as follows: Step 8.1, using the MUSCL reconstruction method to calculate the spatial variation gradient of each variable in the unit grid; obtaining the restricted variable gradient by restricting the variable gradient, and interpolating the variable value at the midpoint of the unit boundary based on the restricted gradient to achieve MUSCL second-order reconstruction; Step 8.2: Obtain the second-order time accuracy by the two-step Runge-Kutta method. The values ​​are: Intermediate variable value for Where i is the number of control unit grids; Ω is the volume of the i-th control unit grid, in m 3 ;q n+1 is the single-width flow at time step n+1, in m 2 / s,q n is the single-width flow at time step n, in m 2 / s; Δt is the time step, in seconds; S(q n ) is the source term at time step n; during the entire solution process, GPU parallel computing technology is introduced to improve the computing efficiency; Step 8.3: Under the spatial step length Δx, the calculation formula for the value of the time step length Δt is as follows: In the formula, Δx min It represents the minimum distance from the center of the grid unit to the corresponding interface, in meters; u is the flow velocity, in meters per second; g is the gravitational acceleration, with a value of 9.81 meters per second 2 ; h is the water depth, in meters.

9. The efficient water eutrophication assessment method of integrated inorganic and organic substances according to claim 8 is characterized in that: The step 9 is specifically as follows: Step 9.1, specific nutritional status calculation formula: In the formula, TLI(C chla ) is the chlorophyll α nutritional status index of phytoplankton; TLI (C TP ) represents the total phosphorus nutritional status index; TLI (C TN ) is the total nitrogen nutritional status index; TLI (C COD ) is the chemical oxygen demand nutrient status index; C chla is the concentration of chlorophyll α, in mg / L; C TP is the total phosphorus concentration, in mg / L; C TN is the total nitrogen concentration, in mg / L; C COD is the chemical oxygen demand concentration, in mg / L; Step 9.2: Using the above calculation formulas for each nutritional status index and considering the weight of each nutritional status index, the formula for the comprehensive nutritional status index is: Where, TLI(∑) is the comprehensive nutrient status index determined by the above environmental factors, reflecting the eutrophication of rivers and lakes; TLI(j) is the nutrient status index of the jth parameter; W j is the corresponding weight value of the jth parameter nutritional status index; Taking the chlorophyll α concentration of phytoplankton as the benchmark, the relevant weight of the jth parameter is obtained, and the normalized formula is: In the formula, is the correlation coefficient between the jth parameter and the chlorophyll α concentration; m is the number of parameters, which is used to evaluate the eutrophication of river and lake water bodies.

Images