CN109472072B

CN109472072B - Seasonal river and underground water interaction prediction method based on river simulation

Info

Publication number: CN109472072B
Application number: CN201811277574.8A
Authority: CN
Inventors: 陆垂裕; 李天辰; 曹国亮; 孙青言; 严聆嘉; 赵勇; 王建华; 张博; 李泽鹏; 秦韬
Original assignee: China Institute of Water Resources and Hydropower Research
Current assignee: China Institute of Water Resources and Hydropower Research
Priority date: 2018-10-30
Filing date: 2018-10-30
Publication date: 2020-09-01
Anticipated expiration: 2038-10-30
Also published as: CN109472072A

Abstract

The invention discloses a seasonal river and underground water interaction prediction method based on river simulation, and an algorithm for automatically determining a simulation path is designed according to characteristics of a river network and is applied to a large-scale complex natural/artificial river network. The river numbering is random, the confluence relation and the diversion relation of a river network do not need to be considered, and the work of a user is facilitated; meanwhile, the invention improves the function of multi-river diversion, and can simulate the quantitative diversion and proportional diversion of a plurality of rivers and the situation that the dividable water volume is more than or less than the total water demand; in addition, the method has higher accuracy in predicting the interaction flow between the seasonal river and the underground water.

Description

Seasonal river and underground water interaction prediction method based on river simulation

Technical Field

The invention belongs to the technical field of automatic processing of hydrological water resource data information, and particularly relates to a design of a seasonal river and underground water interaction prediction method based on river simulation.

Background

The interconversion between surface water and underground water is a hydrological phenomenon ubiquitous in nature, and is one of the basic properties of water resources. In the water resource evaluation and planning management, the close hydraulic connection between surface water and underground water must be paid great attention, and the surface water and the underground water are coupled and considered to provide a solid foundation for the water resource evaluation and planning management. In the interaction system of surface water and underground water, the significance of the research on river-aquifer interaction is the most important, namely the relation between seasonal rivers and underground water. Seasonal rivers mainly refer to intermittent rivers in which the river bed is exposed when the river water is cut off in the dry season and the river water is running in the rich season.

In order to depict the relation between the underground water flow and the seasonal river, the underground river system simulation analysis is needed, at first, the underground water simulation and prediction mostly adopt a simple water balance method and a hydrogeological simulation method, the recognition of boundary conditions usually treats the river as a fixed water head boundary or a given water head boundary, and the situation is not in accordance with the actual situation in many cases. In practical situations, the water level of the river and the flow rate exchanged with the aquifer are dynamically changed, and the complex mechanism is difficult to describe by a traditional method. With the development of computer technology, numerical simulation is applied to underground water analysis and calculation, the calculation theory and the actual measurement of underground water analysis capability are greatly improved, and the dynamic relationship between underground water and surface water is well described. The existing numerical simulation methods of the underground river system mainly comprise a finite difference method, a finite element method, a boundary element method, a finite volume method and the like, and various underground water numerical simulation software including common software such as VisualMODFLOW, FEFlow, GMS and the like continuously emerges for decades based on the simulation methods. At present, a numerical simulation model represented by MODFLOW based on a finite difference method of a unit center solves the problem of partial differential equation definite solution describing various states of underground water by using an iterative solution method.

However, the numerical simulation software based on the finite difference method of the central unit disperses the research area in time and space, substitutes the differential quotient approximation for the derivative quotient on a discrete grid, conditions the differential equation and the definite solution thereof into an algebraic equation, and iteratively solves to obtain the water head value. Since each discrete grid cell contains both input and output terms, the upstream outflow is equal to the downstream inflow, and therefore the grid cells must be ordered in the order of river confluence. Most simulation software requires users to specify the river confluence relation by themselves in this respect: all rivers are sequenced, all rivers in the river system are arranged into a sequence according to the confluence sequence, and the rivers in the sequence are named one by one according to the serial number so as to facilitate the subsequent simulation operation. However, in the complex river network with criss-cross structure, the work is complicated and error-prone, if the river network relationship or the confluence sequence is not properly arranged, the simulation result of the whole program is wrong in calculation, and under the condition that the program is not checked, the generated error is concealed and is difficult to be checked.

In addition, the existing simulation software is not flexible in processing the manual diversion river when a river network is generated, when the river needs to be diverted according to the diversion ratio, a user is required to convert the diversion ratio into a fixed flow value and input the fixed flow value into each diversion item, and the situation that the water diversion amount is smaller than the total water demand of all downstream diversion rivers cannot be processed.

Disclosure of Invention

The invention aims to solve the problems of numerical prediction of interaction between seasonal rivers and underground water by adopting the conventional simulation software, and provides a seasonal river and underground water interaction prediction method based on river simulation.

The technical scheme of the invention is as follows: the method for predicting the interaction between the seasonal river and the underground water based on river simulation comprises the following steps:

and S1, spatially dividing the river network system of the seasonal river into four levels of a river network, a river system, a river and a river reach.

And S2, randomly numbering rivers and river reach.

And S3, acquiring the data information contained in each river and the data information contained in each river reach, and storing the river data information and the river reach data information into different fields of the database.

And S4, constructing a river structure according to the river data information and the river reach data information.

And S5, constructing a river network confluence structure according to all the river structures.

And S6, recognizing the simulation sequence of each river in the river network confluence structure and generating a simulation path.

And S7, calculating the flow between each river in the simulation path and the underground water.

S8, judging whether the flow between the river and the underground water in all stress periods is calculated, if so, entering the step S9, otherwise, returning to the step S4.

And S9, storing the calculation result into a database, and completing the prediction of the interaction between the seasonal river and the underground water.

Further, the river data information in step S3 includes:

the river number, starting with number 1, is stored in the field SEGMID.

The downstream number, representing the downstream river number for that river, is-1, representing no downstream river, stored in field nexttint.

Head-end incoming traffic is stored in field INFLOW.

The diversion status, which is not-1, indicates the river number to be diverted from the current river, and-1 indicates that the current river is not diverted from another river and is stored in the field DIVSEGM.

In the diversion mode, if the current river is diverted from another river, a value of 0 indicates diversion at a predetermined flow rate, and a value of 1 indicates diversion at a diversion ratio, and the value is stored in the field divorpt.

And the water diversion ratio represents the proportion of the river diversion, the value range is 0-1, and the water diversion ratio is stored in a field DIVRATIO.

And calculating the water level, wherein the value of 1 represents that the water level of the river is calculated according to the Manning formula, and the value of 0 represents that the water level is specified by a user and is stored in a field BCALSTAGE.

The river reach data information in step S3 includes:

stress period number, stored in field IPERIOD.

The grid number, which represents the layer, line and column number of the cell in which the current river reach is located, is stored in the fields ILYR, IROW and ICOL, respectively.

The river number is stored in the field SEGMID.

The intra-river number, which indicates the number of the current river reach within the river, is stored in the field CELLID.

Water level, stored in field strstart.

The hydraulic conductivity, which represents the hydraulic conductivity between the current river reach and the underground aquifer, is stored in the field COND.

Riverbed top elevation, stored in field STRTOP.

The riverbed bottom elevation is stored in the field STRBTM.

Riverbed width, stored in field strwd.

The bed slope, stored in field STRSLP.

The Mannich roughness coefficient, stored in field STRNDC.

Further, when the value of the field divorpt is 1, the shunting mode of the river is shunting from the shunted river according to the shunting ratio, and the shunting formula is as follows:

where n is the number of the divided rivers, i represents the ith divided river and i is 1,2_iIndicating the divided flow of the ith divided stream, α_iSetting the diversion ratio for the ith diversion river system, β_iFor the actual split ratio of the ith split stream, QNRIV represents the remaining flow into the downstream stream after splitting.

If the sum of the water diversion proportions of the diversion rivers is less than or equal to 1, the outflow flow is diverted according to the water diversion proportions, and the residual flow is converged into the downstream river; if the sum of the water division ratios is larger than 1, the downstream of the river has no flow, and each divided river is divided according to the water division ratio.

When the value of the field DIVOPT is 0, the shunting mode of the river is to shunt from the shunted river according to the specified flow, and the shunting formula is as follows:

wherein Q_iThe designated flow representing the score of the ith river is positioned in the field INFLOW; if the outflow of the shunted rivers is more than or equal to the specified flow of each shunted river, shunting according to the specified flow, and converging the residual flow into the downstream river; if the outflow of the shunted rivers is less than the designated flow of each shunted river, no flow exists at the downstream of the rivers, and each shunted river shunts according to the designated flow proportion.

Further, the river structure constructed in step S4 includes:

and the upstream unit pointer array is used for storing pointers of all upstream river units converged to the local river unit.

And the diversion unit pointer array is used for storing pointers of the river units to be diverted by the river units.

And the downstream unit pointer is used for storing the pointer of the downstream river flow unit of the river flow unit.

The river cell data array comprises river cell number data, river cell attribute data and river cell diversion data.

And the river reach data array is used for storing the data information of all river reaches in the river unit.

Further, step S5 includes the following substeps:

s51, selecting the first river structure according to the serial number of the river.

S52, judging whether the river structure has a downstream river structure, if so, going to step S53, otherwise, going to step S54.

S53, the downstream cell pointer of the river structure is pointed to the upstream cell pointer array of the downstream river structure, and the process advances to step S55.

S54, marking the river structure as a top-level river outlet, adding the top-level river outlet into the top-level river outlet array, and proceeding to the step S55.

S55, judging whether all the river structure bodies have completed the pointer relation setting, if yes, entering step S56, otherwise, selecting the next river structure body, and returning to step S52.

S56, selecting the first river structure according to the serial number of the river.

And S57, judging whether the river structure is a diversion river structure, if so, entering the step S58, and otherwise, entering the step S59.

S58, the pointer in the upstream cell pointer array of the divided stream structure is pointed to the divided stream cell pointer array of the upstream stream structure, and the process advances to step S59.

S59, judging whether all the river structural bodies have finished the setting of the stream dividing relationship, if so, finishing the construction of the river network confluence structure, and going to step S6, otherwise, selecting the next river structural body and returning to step S57.

Further, the process of identifying the simulation order of each river in step S6 includes river scanning and river traversal.

River scanning comprises the following steps:

a1, initial set scan number isuhmscan 0.

A2, initializing the scan termination flag ifiish ═ 0.

A3, setting the river system number i in the river network confluence structure to 1.

A4, calling a recursive function RoutetTest (i) to traverse the ith river system, and acquiring a return value of the function RoutetTest (i); if the return value of the function routtest (i) is 0, the simulation sequence of each river in the current river system is determined, and if the return value of the function routtest (i) is-1, the simulation sequence of a river in the current river system cannot be determined in the current river system scanning.

A5, adding the return value of the current lineage function RoutetTest (i) to the scanning termination identifier iFinish.

A6, adding 1 to the river series number i.

A7, judging whether i is larger than the number of river systems in the river network confluence structure, if yes, entering the step A8, otherwise, returning to the step A4.

A8, judging whether the scanning termination identifier iFinish is equal to 0, if so, indicating that the simulation sequence of all rivers in the river system is determined, ending the river system scanning process, otherwise, entering the step A9.

A9, add 1 to the number of scans isuhmscan.

A10, judging whether the scanning times iSumScan are larger than the total number of river units in the river network confluence structure, if so, indicating that the physical structure of the river network is unreasonable, simulating the abnormality, ending the river system scanning process, otherwise, returning to the step A2.

River system traversal comprises the following steps:

and B1, setting the ordinal number of the currently traversed river structure as j.

B2, judging whether the simulation sequence of the j river structure is determined in the previous river system scanning process, if so, setting the return value of a recursive function RoutetTest (j) to be 0, and entering the step B13, otherwise, entering the step B3.

B3, judging whether the upstream cell pointer array of the jth river structure has a value, if so, indicating that a branch river structure exists upstream of the jth river structure, and entering step B4, otherwise, indicating that no branch river structure exists upstream of the jth river structure, and entering step B11.

B4, initializing the setting result value iResult ═ 0.

B5, initializing the number k of upstream branch river structure to be 1.

B6, calling a recursive function routtest (k) to traverse the kth upstream tributary river structure.

B7, acquiring the return value of the function RoutetTest (k) and accumulating the return value into a result value iResult; if the return value of the function routtest (k) is 0, it indicates that the simulation order of the upstream branch river structure is determined, and if the return value of the function routtest (k) is-1, it indicates that the simulation order of the upstream branch river structure is not determined.

B8, adding 1 to the structure ordinal number k of the upstream branch river.

B9, judging whether k is larger than the number of the upstream branch river structure bodies, if so, entering a step B10, otherwise, returning to the step B6.

B10, judging whether the result value iResult is equal to 0, if so, indicating that the simulation sequence of the jth river structure is determined, setting the return value of routtest (j) to 0, and entering step B13, otherwise, indicating that the simulation sequence of the jth river structure cannot be determined, setting the return value of routtest (j) to-1, and entering step B13.

B11, determining whether the diversion unit pointer array of the jth river structure has a value, if so, indicating that the jth river structure is diverted from other river structures, and going to step B12, otherwise, indicating that the jth river structure is not diverted from other river structures, and is the upstream first-level river of the current river system, further determining the simulation sequence of the jth river structure, setting the return value of routtest (j) to 0, and going to step B13.

B12, judging whether the simulation sequence of the j river structure shunted by the j river structure is determined in the previous river system traversal process, if so, determining the simulation sequence of the j river structure, setting the return value of RoutetTest (j) to be 0, and entering the step B13; otherwise, the simulation sequence of the jth river structure cannot be determined, the return value of routtest (j) is set to-1, and the process proceeds to step B13.

B13, outputting the return value of routtest (j), and ending the traversal of the jth river structure.

Further, the flow rate between the river and the groundwater in step S7 is calculated by the formula:

wherein QRIV' represents a flow rate between the river and the groundwater and is stored in the field STREAM, and a positive value thereof represents a leakage amount of the river to the groundwater aquifer and a negative value thereof represents a discharge amount of the groundwater aquifer to the river; CRIV represents the hydraulic conductivity coefficient between river-underground aquifer, stored in field COND; h is_i,j,kThe calculated water head of a unit where the river reach is located is represented, RBOT represents the elevation of a certain point at the bottom of a river bed, and the elevation is stored in a field STRBTM; HRIV indicates the water level of the river reach and is stored in the field STRSTAGE.

Further, the method for determining the water level HRIV of the river reach comprises the following steps:

if the value of the field BCALSTAGE is 0, the water level HRIV of the river reach is specified by a user; if the value of the field BCALSTAGE is 1, the river reach water level is calculated according to a Manning formula, wherein the calculation formula is as follows:

wherein Q is the inflow rate of each river section of the river, n is a Manning roughness coefficient and is stored in a field STRNDC, and c is the hydraulic conductivity between the river section and the underground aquifer and is stored in a field COND; w is the bed width of the river reach and is stored in field strwd; s is the bed slope of the river reach and is stored in field STRSLP.

The invention has the beneficial effects that:

(1) the river numbering is carried out randomly in the invention, the confluence relation and the diversion relation of the river network do not need to be considered, the work of a user is facilitated, after the user numbers all the rivers in a research area, only the number of the downstream river channel into which each river flows needs to be further pointed out, and then the rivers in the whole river network can be divided according to the confluence sequence and named one by one according to the sequence by an algorithm for automatically searching the simulation sequence of each river in the river network, so that the requirement of quantitative analysis is met.

(2) The invention improves the function of multi-river diversion, and can simulate the quantitative diversion and proportional diversion of a plurality of rivers and the situation that the dividable water volume is more than or less than the total water demand.

(3) The method has higher accuracy in predicting the interaction flow between the seasonal river and the underground water, and the relative error is only 0.005 percent compared with the actual situation.

Drawings

Fig. 1 is a flowchart illustrating a seasonal river and groundwater interaction prediction method based on river simulation according to an embodiment of the present invention.

Fig. 2 is a schematic diagram illustrating a relationship between the groundwater differential unit and a river or a river reach according to an embodiment of the present invention.

Fig. 3 is a schematic view of the interior of a river structure according to an embodiment of the present invention.

Fig. 4 is a schematic view of a river network confluence structure according to an embodiment of the present invention.

Fig. 5 is a flow chart illustrating a structure for constructing a river network confluence structure according to an embodiment of the present invention.

Fig. 6 is a flow chart of river scanning according to an embodiment of the present invention.

Fig. 7 is a flow chart of river traversal provided by the embodiment of the present invention.

FIG. 8 is a schematic view illustrating the underground water level dropping below the bottom layer of the riverbed according to the embodiment of the present invention.

Fig. 9 is a schematic view of the bottom surface of the riverbed bottom layer with the groundwater level higher than the riverbed bottom layer according to the embodiment of the invention.

Fig. 10 is a schematic view of a river network structure and a diversion situation provided in the embodiment of the present invention.

Fig. 11 is a schematic view illustrating a river network confluence relationship according to an embodiment of the present invention.

Fig. 12 is a schematic view illustrating an outflow rate of a river network unit according to an embodiment of the present invention.

Fig. 13 is a schematic view illustrating leakage of a river network unit according to an embodiment of the present invention.

Detailed Description

Exemplary embodiments of the present invention will now be described in detail with reference to the accompanying drawings. It is to be understood that the embodiments shown and described in the drawings are merely exemplary and are intended to illustrate the principles and spirit of the invention, not to limit the scope of the invention.

The embodiment of the invention provides a method for predicting interaction between a seasonal river and underground water based on river simulation, which comprises the following steps S1-S9 as shown in figure 1:

The river network refers to a surface water channel system formed by all cross rivers in the research area. The river in the river network may flow out through a plurality of watershed outlets, thus dividing the river network into a plurality of river systems, each river system corresponding to one watershed outlet. The river system is composed of a plurality of rivers, wherein the rivers are a certain channel through which surface water flows, and the river system has the concept of upstream and downstream, water flows to downstream through the upstream of each river and is merged into the river channel at the next stage, so that the inflow of the upstream of the river at the next stage is formed, and the water flows out from the downstream outlet (also called a watershed outlet) of the lowest river of the river system through the step-by-step confluence of the rivers.

In the embodiment of the invention, each river is a surface water channel unit which is independent from each other, and in order to simulate the mutual coupling action between the surface water and the underground water of the river, the river needs to be divided into river sections according to the space distribution relation between the underground water differential grids of the river. The river reach is a section of a river distributed in a grid unit of underground water, and is a basic unit for simulating the interaction of surface water and underground water of the river. The water quantity in the river flows from the upstream to the downstream, so that the river is divided into sections from the upstream to the downstream according to the underground water differential grid units, and the river sections are numbered in sequence, and the relative sequence of the river section numbers from small to large represents the flowing direction of the water flow. Fig. 2 shows a schematic relationship among rivers, groundwater differential grid cells, and river segments, and fig. 2 shows 7 rivers, in which the 1 st river is divided into 2 segments and the 2 nd river is divided into 4 segments.

And S2, randomly numbering rivers and river reach.

Each river in the study area has a unique number for distinguishing different rivers, and the numbers of the rivers are numbered continuously from 1, so that the maximum river number is the number of rivers in the study area. In the embodiment of the invention, the number of the river has no other meaning except for identifying the river, so that the numbers of all rivers in the research area can be randomly disordered and have high freedom degree of numbering. Other prior art methods usually require the number of the river and also indicate the front-back sequence of the river in the simulation process, such as MODFLOW and the like, which means that a user must identify the confluence relation and the diversion relation between the rivers in the river network and then sequentially number the rivers, and the sequence of the numbers represents the front-back sequence of the rivers in the simulation process, which is not only a heavy workload for a more complex river network, but also is easy to make mistakes. According to the invention, when numbering rivers, a user does not need to consider the confluence relation and the diversion relation of a river network, so that the work of the user is facilitated. After the user numbers all rivers in the research area, the user only needs to further indicate the number of the downstream river channel into which each river flows, and in the subsequent steps, a river network confluence structure can be automatically established through a corresponding algorithm, and the simulation sequence of each river in the river network is identified. The numbering of the 7 rivers in the study area is shown in fig. 2.

Each river is divided into different river segments according to the distribution of each river in the underground water differential grid, and the division and numbering rules of the river segments are shown in fig. 2. Cell (1,3) represents a cell located at row 1, column 3. River 1 has two sections, starting at cell (1,3) and ending at cell (2,3) along the direction of river flow. Some of which is introduced into river 2 and the remainder of which enters river 3 from river 1, and river 3 has 4 segments. River 2 and river 4 meet in cell (5,3) to form river 5 containing 2 river segments, and river 3, river 5 and river 6 meet in cell (5,4) to form river 7. Two small portions of rivers in stream 3 and stream 6 are not included in the numbering plan because they only pass through the corners of cells (2,4) and (3, 5). In addition, multiple segments of different rivers may be assigned to the same model cell, such as the second segment of the 5 th river; the 4 th stretch of the 3 rd river; the 5 th stretch of the 6 th river; the 1 st river segment of the 7 th river is located in the unit (5, 4).

The river includes data information as shown in table 1, and different data information is stored in different fields of the database.

TABLE 1

In the embodiment of the invention, when the value of the field DIVOPT is 1, the shunting mode of the river is to shunt from the shunted river according to the water diversion ratio, and the shunting formula is as follows:

wherein n is the number of the shunted rivers (n is more than or equal to 1), i represents the ith shunted river and i is 1,2_iIndicating the divided flow of the ith divided stream, α_iSetting the diversion ratio for the ith diversion river system, β_iFor the actual split ratio of the ith stream, QNRIV indicates the remaining flow into the downstream stream after splitting in an embodiment of the present invention, the split ratio α for each split stream is set by the system_iThe sum may be greater than 1, thus introducing an actual split ratio β_iAs the actual water diversion ratio of the ith diversion river in the simulation.

The river reach contains data information as shown in table 2, and the information data is the basis of simulating the hydraulic connection between surface water and underground water of the river.

TABLE 2

On the basis of the preparation of river and river reach data information, a data structure body of a river unit can be designed, so that the organization and the use of information data are convenient in the simulation process, the internal structure of the river structure body is shown in fig. 3, and each river structure body comprises:

On the basis of the river structure, a confluence structure of the whole river network of the research area can be constructed by using the target pointer of the river structure. Generating a river structure object according to the data information of each river and the river reach thereof, and connecting different river structures by using the upstream and downstream relations of the river, the river diversion relation and the river structure object pointer to form a river network convergence structure in the simulation, as shown in fig. 4.

In the process of constructing the river network confluence structure, two factors of nature and manpower need to be considered, namely the upstream and downstream relations of different rivers in the river network, namely the topological structure of the natural river network; and the other is the manual diversion relation among different rivers, namely the water transfer relation under human intervention. Under natural conditions without human intervention, the river is always converged by the upstream primary river to the downstream basin outlet. As shown in fig. 5, the construction process of the river network confluence structure includes the following steps S51 to S59:

Before step S51, it is first necessary to extract river information and attribute information from all the fields in the parsing database about seasonal river boundary data, and determine a river number. After the serial number continuity is checked, the serial number information of the river is stored by adopting a hash table structure, so that the inquiry is convenient during subsequent analysis and calculation.

The specific method for judging whether a certain river structure is a divided river structure is as follows: when the head inflow value in the river structure unit data is not 0 and the diversion river number is not-1, the river structure is the diversion river structure. The head-end INFLOW and the shunt river number are respectively obtained by analyzing the INFLOW and DIVSEGM fields in the database.

The automatic identification of the simulation sequence of each river in the river network confluence structure is the core of the invention. In the embodiment of the invention, firstly, a simulation sequence array is set, and the object pointer of the river structure body is arranged at the tail of the simulation sequence array when the simulation sequence of a river is identified. Thus, after the simulation sequences of all rivers are identified, the sequence of each element in the simulation sequence array represents the simulation sequence of each river.

The algorithm flow for identifying the river simulation sequence includes two major steps: river scanning and river traversal. If the rivers in the river network are not connected by artificial water diversion, the simulation sequence of each river can be directly determined by traversing each river system from upstream to downstream; if artificial water diversion connection exists among rivers in the river network, a simulation path needs to be determined through a method combining river system traversal and river system scanning.

Since the river network confluence structure has a complicated water diversion relationship, if a river to be diverted is not simulated before the river simulation, the flow rate of the diverted river cannot be determined during the simulation, and the diversion process cannot be realized. Because the diversion relation is completely independent of the natural step-by-step convergence relation, no rule can be stated, and the final reasonable river simulation sequence can be determined only by continuously scanning and testing each river system in the river network.

And calling the recursive function RoutetTest one by one for all the river systems in the river network convergence structure so as to complete the traversal of all the river systems, which is called as one-time scanning. When no inter-river diversion relationship exists in the river network, the confluence sequence of each river system is independent, so that only 1 scanning is needed, and the simulation sequence of each river in the river network system can be determined. However, when there is a water diversion relationship between rivers, a river in a certain river system may need to wait for a river in another river system to be simulated before determining the upstream water amount, so that a scanning process of the river system cannot complete the sequencing of all river simulation sequences in the river system. In the embodiment of the present invention, an algorithm for repeatedly scanning a river network is provided to ensure that the simulation sequence of each river in each river system can be reasonably determined, as shown in fig. 6, the river system scanning includes the following steps a 1-a 10:

a1, initial set scan number isuhmscan 0.

A2, initializing the scan termination flag ifiish ═ 0.

A6, adding 1 to the river series number i.

A9, add 1 to the number of scans isuhmscan.

Each time the scan number isuhmscan is updated, it needs to be checked. Since the total number of river cells in the river network is fixed and the simulation order of at least one more river can be determined per scan, there is a threshold for the number of scans isuhmscan. If the value of the scanning times iSUMScan is larger than the total number of the river units in the river network, the simulation is abnormal, the physical structure of the river network is unreasonable, the simulation is quitted, and the river scanning process is ended.

Each river system in the river network has a top river outlet (the most downstream river unit), the river system traversing process adopts a recursive calling method, and the upstream river structure is traversed from the river basin outlet of each river system until the most upstream river flow (the river which is not shunted by other river structures) is found, and then the river system traversing process is sequenced from the primary river to the downstream. As shown in FIG. 7, river traversal comprises the following steps B1-B13:

B4, initializing the setting result value iResult ═ 0.

B5, initializing the number k of upstream branch river structure to be 1.

B8, adding 1 to the structure ordinal number k of the upstream branch river.

The flow calculation formula between the river and the underground water is as follows:

wherein QRIV' represents a flow rate between the river and the groundwater and is stored in the field STREAM, and a positive value thereof represents a leakage amount of the river to the groundwater aquifer and a negative value thereof represents a discharge amount of the groundwater aquifer to the river; CRIV represents river-underground water contentHydraulic conductivity coefficients between the layers, stored in field COND; h is_i,j,kThe calculated water head of a unit where the river reach is located is represented, RBOT represents the elevation of a certain point at the bottom of a river bed, and the elevation is stored in a field STRBTM; HRIV indicates the water level of the river reach and is stored in the field STRSTAGE.

Calculating the flow QRIV' requires the use of three parameters: the water level HRIV of the river reach, the hydraulic conductivity coefficient CRIV between the river and the underground aquifer and the elevation position RBOT at the bottom of the river bed, and the parameters in the QRIV' expression (3) need to be split and substituted into a difference equation based on water quantity equalization.

When h is generated_i,j,kWhen RBOT is not more than equal to, h is shown in figure 8_i,j,kHas fallen below the bottom surface of the bed bottom layer, forming an unsaturated zone thereunder. If the head of water at a point at the bed bottom is simply considered as the elevation of that point, i.e., RBOT, assuming that the bed layer itself remains saturated, the flow between the river and the groundwater is a constant value. the-CRIV (HRIV-RBOT) needs to be added to the differential term RHS at the time of calculation; when h is generated_i,j,k>RBOT, as shown in FIG. 9, h_i,j,kThe flow rate between the river and the underground water is proportional to the head difference between the river and the underground aquifer and is calculated by the head h_i,j,kAs variables. During calculation, the-CRIV is added into the difference term HCOF, and the-CRIV × HRIV is added into the difference term RHS, so that the QRIV' split is substituted into the difference equation based on water balance.

And (4) connecting the difference equations obtained by all the calculation units, and solving by adopting an iterative method. In the iteration process, the result of each iteration is processed and used for the next calculation. Different algorithms have different processing methods, and under normal conditions, the water head change after each iteration is gradually reduced, and finally convergence is achieved. This completes the calculation of the head for a period of time. Convergence is usually determined by a predefined convergence index, and is referred to as convergence when the maximum head difference calculated in two iterations is less than the convergence index. The water head value at the end of each time period is obtained from the initial water head step by step, and is used as the initial value of the next time period, and the process is repeated until the required time is ended, so that the flow between each river and the underground water in the simulation path can be obtained.

The method for determining the water level HRIV of the river reach comprises the following steps:

In the embodiment of the invention, no matter the designated water level calculation or the automatic water level calculation is carried out according to the Manning formula, the upstream inflow of one river reach is subtracted by the leakage amount of the current river reach (or added with the leakage amount of underground water to the current river reach), and then the downstream outflow serving as the current river reach flows to the next river reach to become the upstream inflow of the next river reach. This in turn allows the river discharge to be calculated from the most upstream section of the river to the most downstream section of the river, which in turn will become the upstream inflow of the next river, and this process will continue until the river discharge leaves the river network system.

The effect of the present invention is further illustrated by the following specific example:

as shown in fig. 10, there were 9 rivers in the river network, and the river numbers were not in the order of confluence of the river network, and the respective river manning roughness coefficients were all 0.013.

Wherein, the river 1 has 6 river reach with the length of 2333.21m, 5203.9m, 3093.2m, 5303.9m, 5589.5m and 2748.8 m; the elevation of the top of the riverbed is-0.5 m; the elevation of the bottom of the riverbed is-1.1 m; the width of the riverbed is 12 m; the slope of the river bed is 0.0012.

The river 2 has 7 river reach with lengths of 4302.4m, 5003.7m, 5003.7m, 5003.7m, 5003.7m, 5003.7m and 3279.8 m; the elevation of the top of the riverbed is-2 m; the elevation of the bottom of the riverbed is all-3 m; the width of the riverbed is 20 m; the slope of the riverbed is 0.001.

The river 3 has 10 river sections with the lengths of 44615.4m, 6349.3m, 6350.4m, 6350.4m, 6350.4m, 6350.4m, 6350.4m, 6350.4m, 6350.4m and 3657.2 m; the elevation of the top of the riverbed is-2.5 m; the elevation of the bottom of the riverbed is-3.3 m; the width of the riverbed is 20 m; the slope of the riverbed is 0.001.

The river 4 has 4 river reach with lengths of 5584.0m, 2750.6m, 4756.1m and 6507.3 m; the elevations of the top of the riverbed are all-2.6 m; the elevation of the bottom of the riverbed is-3.4 m; the width of the riverbed is 20 m; the slope of the riverbed is 0.001.

The river 5 has 9 river sections with the lengths of 2872.5m, 4636.0m, 3969.2m, 5190.0m, 5190.0m, 3572.2m, 4983.6m, 5190.0m and 5190.0 m; the elevation of the top of the riverbed is-3 m; the elevation of the bottom of the riverbed is-4 m; the width of the riverbed is 20 m; the slope of the riverbed is 0.001.

The river 6 has 4 river reach with lengths of 2647.4m, 6812.4m, 8235.4m and 4735.2 m; the elevations of the top of the riverbed are all-2.3 m; the elevation of the bottom of the riverbed is-3.1 m; the width of the riverbed is 20 m; the slope of the riverbed is 0.001.

River 7 has 8 river reach with lengths of 5463.0m, 8445.2m, 8445.2m, 8445.2m, 8445.2m, 3686.0m, 8445.2m and 4759.2 m; the elevation of the top of the riverbed is-3.2 m; the elevation of the bottom of the riverbed is-3.8 m; the width of the riverbed is 20 m; the slope of the riverbed is 0.001.

The river 8 has 4 river reach with lengths of 7161.3m, 13369.2m, 16279.3m and 6149.1 m; the elevation of the top of the riverbed is-5 m; the elevation of the bottom of the riverbed is-9 m; the width of the riverbed is 20 m; the slope of the riverbed is 0.001.

River 9 has 10 river reach lengths of 2781.2m, 3001.0m, 2321.0m, 2543.0m, 2442.0m, 2797.0m, 2611.0m, 2351.0m, 2124.0m and 2100.0 m; the elevation of the top of the riverbed is-1 m; the elevation of the bottom of the riverbed is-2 m; the width of the riverbed is 12 m; the slope of the river bed is 0.0012.

The permeability coefficient of the unit grid along the row direction is 8 m/s; the vertical hydraulic conductivity is 0.00002; the initial head of the grid cell is-1 m. Initial flow rate of river 1 was 12000m³Min, initial flow of river 8 is 10000m³The water flow rate of 7 parts of river 3 from river is 4000 m/min³And/min. River 5, river 9 and river 2 the proportion of water divided from river 1 is 3:6:1 respectively. The space unit of the simulation is m; the time unit is day.

Firstly, steps S2-S5 are adopted to construct a river network confluence structure, and command keywords in a database are searched, wherein the command keywords comprise SEGMID, NEXTUNIT, DIVSEGM and other fields. As shown in table 3, the downstream river of river 1 is river 2; the downstream river of river 2 is river 4; the downstream river of river 4 is numbered-1, i.e., river 4 has no downstream river, being the top level river outlet. Thus, each river structure is analyzed, and a river network confluence structure is established.

TABLE 3

Setting a water regulation relation on the basis of a natural river network: stream 3 is split from stream 7; river 5 and river 9 are branched from river 1; the river no-diversion corresponding to the DIVSEGM field of-1. All the connected river units are found and pointer relations are set, and the formed river network confluence relations are shown in fig. 11.

The simulated path is then automatically searched through the river traversal and river scanning in step S6. The river network confluence structure has two river systems in common, namely a left river system (comprising

rivers

1,2, 3 and 4) and a right river system (comprising

rivers

5, 6, 7, 8 and 9). Traversing from the top level river basin outlet of the river system to the upstream river.

And performing first river system scanning, and traversing the left river system. From the river outlet, the river 4 is traversed, the

upstream branch rivers

2 and 3 are found, and the river 2 is recursively called first. And traversing the river 2, finding the upstream river 1, and performing recursive calling on the river 1. And traversing river 1, wherein the river 1 is the initial river upstream of the left river system and has no upstream river, adding the river into the simulation sequence group and returning to the value 0, ending the recursive call of the river 1 and returning to the recursive call function of the river 2. River 2 gets a return value of 0, upstream of which the simulation order has been determined, adds river 2 to the set of simulation orders and returns a value of 0. The recursive call of river 2 is ended and the recursive call function of river 4 is returned. River 4 takes the first upstream tributary return value 0 and then traverses the second upstream river 3. River 3 is diverted from river 7 and river 7 has no determined simulation order, the left river system traversal is stopped and returns to-1. The second upstream return value of river 4 is-1, the sum of the two return values is-1, and the return value is returned to the left river system. Then, the right river system is traversed. From the river outlet, the river 6 is traversed, the

upstream branch rivers

9 and 7 are found, and the river 9 is recursively called first. Traversing river 9, river 9 is diverted from river 1 and river 1 has determined the simulation order, adding river 9 to the set of simulation sequences and returning a value of 0. The recursive call of river 9 is ended and the recursive call function of river 6 is returned. River 6 gets a first upstream return value of 0 and then traverses a second upstream river 7. When the river 7 is traversed and

upstream branch rivers

5 and 8 are found, the river 5 is recursively called first. Traversing river 5, river 5 is shunted from river 1 and river 1 has determined the simulation order, adding river 5 to the set of simulation sequences and returning a value of 0. The recursive call of river 5 is ended and the recursive call function of river 7 is returned. River 7 gets the first upstream return value of 0 and then traverses upstream river 8. Traversing river 8, river 8 being the upstream top river, adds it to the set of simulated sequences and returns a value of 0. River 7 obtains a second upstream return value of 0, with the sum of all upstream returns being 0, i.e., all upstream rivers have determined the simulation order, adds river 7 to the set of simulation orders and returns a value of 0. The recursive call of river 7 is ended and the recursive call function of river 6 is returned. River 6 gets a second upstream return value of 0, all of which have determined the simulation order, adds river 6 to the set of simulation sequences and returns 0 to the right river system, at which point all of the rivers of the right river system have determined the simulation order, completing the first river system scan.

The sum of the two river return values is-1, and a second river scan is performed. And traversing the left river system, traversing the river 4, determining the simulation sequence of the upstream river 2 to return a value 0, ending the recursive call of the river 2, and returning to the recursive call function of the river 4. Then traverse a second upstream river 3, river 3 is shunted from river 7 and river 7 has determined the simulation order, adding river 3 to the set of simulation sequences and returning a value of 0. River 4 gets a second upstream return of 0, with the sum of all upstream returns being 0, adds river 4 to the set of simulated sequences and returns 0 to the left river system. If the return value of the right river is 0, the right river is not traversed, and the second river scanning is completed.

The sum of the return values of the two river systems is 0, the simulation sequence of all the rivers in the river network is determined, and the scanning is stopped.

Through the steps of S7-S9, the flow between the river and the underground water is predicted, and the prediction result is analyzed, so that the technical effects of the invention are as follows:

(1) the capability of automatically searching for the simulation path is embodied. Through the analysis of the river network, the rivers in the whole river network are numbered one by one according to the confluence sequence, and the result is shown in table 4, wherein the IORDER field is an analog sequence, and the UNITID field is the original serial number of the river. Therefore, the original numbering of the rivers is random, the confluence relation and the diversion relation of the river network do not need to be considered, and the work of users is facilitated.

TABLE 4

IPERIOD	IORDER	UNITID
				0	1	1
0	2	2
			0	3	9
0	4	5
			0	5	8
0	6	7
			0	7	3
0	8	6
			0	9	4

(2) The function of quantitatively shunting and proportionally shunting a plurality of rivers is realized. River 5, river 9, and river 2 the flow rates of the three rivers divided from river 1 are respectively: 1979.95m³/min、3959.89m³/min、659.99m³And/min. The ratio is 3:6:1 respectively; the discharge rate of the river 7 is 12916.03m³Min, 4000m through river 3³8915.73m after flow rate of/min³The/min water flows into the downstream river 6.

(3) The accuracy of the simulation prediction result is fully demonstrated. The outflow and leakage of each unit grid of the river network are shown in fig. 12 and 13, and the absolute water balance error of the whole area is-0.822 m³(ii) a The error of the relative water balance is 0.005%. Wherein the leakage amount of the seasonal river to the underground water system is 14303.2m³(ii) a The discharge amount of the underground water system to the seasonal river is 27.9m³。

It will be appreciated by those of ordinary skill in the art that the embodiments described herein are intended to assist the reader in understanding the principles of the invention and are to be construed as being without limitation to such specifically recited embodiments and examples. Those skilled in the art can make various other specific changes and combinations based on the teachings of the present invention without departing from the spirit of the invention, and these changes and combinations are within the scope of the invention.

Claims

1. The method for predicting the interaction between the seasonal river and the underground water based on river simulation is characterized by comprising the following steps of:

s1, dividing a river network system of the seasonal river into four levels of a river network, a river system, a river and a river reach on the space;

s2, randomly numbering rivers and river reach;

s3, acquiring data information contained in each river and data information contained in each river reach, and storing the river data information and the river reach data information into different fields of a database;

s4, constructing a river structure according to the river data information and the river reach data information;

s5, constructing a river network confluence structure according to all the river structure bodies;

s6, identifying the simulation sequence of each river in the river network confluence structure to generate a simulation path;

s7, calculating the flow between each river and the underground water in the simulation path;

s8, judging whether the flow between the rivers and the underground water in all stress periods is calculated, if so, entering a step S9, otherwise, returning to the step S4;

s9, storing the calculation result into a database to complete the prediction of the interaction between the seasonal river and the underground water;

the river data information in step S3 includes:

the river number, starting from number 1, is stored in field SEGMID;

a downstream number, representing the downstream river number of the river, the value-1 representing no downstream river, stored in the field nexttint;

head end incoming traffic, stored in field INFLOW;

the shunting condition, if the value is not-1, the current number to be shunted by the current river is represented, and if the value is-1, the current river is not shunted by other rivers and is stored in a field DIVSEGM;

if the river is shunted from other rivers, the value is 0 to indicate shunting according to the specified flow, and 1 to indicate shunting according to the water diversion ratio, and the shunting is stored in a field DIVOPT;

the water diversion ratio represents the proportion of the river diversion, the value range is 0-1, and the water diversion ratio is stored in a field DIVRATIO;

calculating the water level, wherein the value of 1 represents that the water level of the river is calculated according to a Manning formula, and the value of 0 represents that the water level is specified by a user and is stored in a field BCALSTAGE;

the river reach data information in step S3 includes:

a stress period number stored in the field IPERIOD;

grid numbers, which represent the layer, line and column numbers of the cell in which the current river reach is, and are respectively stored in fields ILYR, IROW and ICOL;

the serial number of the belonged river is stored in a field SEGMID;

the number in the river, which represents the number of the river reach in the river, is stored in a field CELLID;

water level, stored in field STRSTAGE;

the hydraulic conductivity, which represents the hydraulic conductivity between the current river reach and the underground aquifer, is stored in the field COND;

riverbed top elevation, stored in field STRTOP;

riverbed bottom elevation, stored in the field STRBTM;

riverbed width, stored in field strwd;

a bed slope, stored in field STRSLP;

a Mannich roughness coefficient stored in field STRNDC;

when the value of the field DIVOPT is 1, the shunting mode of the river is to shunt from the shunted river according to the water diversion ratio, and the shunting formula is as follows:

where n is the number of the divided rivers, i represents the ith divided river and i is 1,2_iIndicating the divided flow of the ith divided stream, α_iSetting the diversion ratio for the ith diversion river system, β_iFor the actual split ratio of the ith split stream, QNRIV represents the remaining flow into the downstream stream after splitting;

if the sum of the water diversion proportions of the diversion rivers is less than or equal to 1, the outflow flow is diverted according to the water diversion proportions, and the residual flow is converged into the downstream river; if the sum of the water diversion proportions is larger than 1, no flow exists at the downstream of the river, and each diversion river is shunted according to the water diversion proportions;

2. The method for predicting seasonal river interaction with groundwater according to claim 1, wherein the river structure constructed in the step S4 includes:

the upstream unit pointer array is used for storing pointers of all upstream river units converged to the local river unit;

the diversion unit pointer array is used for storing pointers of the river units to be diverted by the river units;

the downstream unit pointer is used for storing the pointer of the downstream river flow unit of the river unit;

the river cell data array comprises river cell number data, river cell attribute data and river cell diversion data;

3. The method for predicting seasonal river interaction with groundwater according to claim 2, wherein the step S5 includes the following substeps:

s51, selecting a first river structure according to the serial number of the river;

s52, judging whether the river structure has a downstream river structure, if so, entering a step S53, otherwise, entering a step S54;

s53, pointing the downstream cell pointer of the river structure to the upstream cell pointer array of the downstream river structure, and going to step S55;

s54, marking the river structure as a top-level river outlet, adding the top-level river outlet into a top-level river outlet array, and entering the step S55;

s55, judging whether all the river structural bodies have finished the setting of the pointer relation, if so, entering a step S56, otherwise, selecting the next river structural body, and returning to the step S52;

s56, selecting a first river structure according to the serial number of the river;

s57, judging whether the river structure is a diversion river structure, if so, entering a step S58, otherwise, entering a step S59;

s58, enabling the pointer in the pointer array of the upstream unit of the diversion river structure to point to the pointer array of the diversion unit of the upstream river structure, and going to the step S59;

4. The method for predicting seasonal river and groundwater interaction according to claim 3, wherein the step S6 of identifying the simulation sequence of each river includes river scanning and river traversal;

the river scanning comprises the following steps:

a1, setting the initial scanning times isuhmscan to 0;

a2, setting the scanning termination identifier ifinishe equal to 0;

a3, setting the number i of river series in the river network converging structure as 1;

a4, calling a recursive function RoutetTest (i) to traverse the ith river system, and acquiring a return value of the function RoutetTest (i); if the return value of the function routtest (i) is 0, the simulation sequence of each river in the current river system is determined, and if the return value of the function routtest (i) is-1, the simulation sequence of a river in the current river system cannot be determined in the current river system scanning;

a5, accumulating the return value of the current river system function RoutetTest (i) into a scanning termination identifier iFinish;

a6, adding 1 to the river series number i;

a7, judging whether i is larger than the number of river systems in the river network confluence structure, if so, entering the step A8, otherwise, returning to the step A4;

a8, judging whether the scanning termination identifier iFinish is equal to 0, if so, indicating that the simulation sequence of all rivers in the river system is determined, ending the river system scanning process, otherwise, entering the step A9;

a9, adding 1 to the scanning times iSUMScan;

a10, judging whether the scanning times iSumScan are larger than the total number of river units in the river network confluence structure, if so, indicating that the physical structure of the river network is unreasonable, simulating the abnormality, ending the river system scanning process, otherwise, returning to the step A2;

the river system traversal comprises the following steps:

b1, setting the ordinal number of the currently traversed river structure as j;

b2, judging whether the simulation sequence of the jth river structure is determined in the previous river scanning process, if so, setting the return value of a recursive function RoutetTest (j) to be 0, and entering the step B13, otherwise, entering the step B3;

b3, judging whether the upstream cell pointer array of the jth river structure has a value, if so, indicating that a branch river structure exists upstream of the jth river structure, and entering a step B4, otherwise, indicating that no branch river structure exists upstream of the jth river structure, and entering a step B11;

b4, setting the result value iResult equal to 0;

b5, initializing the number k of upstream branch river structure as 1;

b6, calling a recursive function RoutetTest (k) to traverse the kth upstream tributary river structure;

b7, acquiring the return value of the function RoutetTest (k) and accumulating the return value into a result value iResult; if the return value of the function routtest (k) is 0, the simulation sequence of the upstream tributary river structure is determined, and if the return value of the function routtest (k) is-1, the simulation sequence of the upstream tributary river structure cannot be determined;

b8, adding 1 to the sequence number k of the upstream branch river structure;

b9, judging whether k is larger than the number of the upstream branch river structure bodies, if so, entering a step B10, otherwise, returning to the step B6;

b10, judging whether the result value iResult is equal to 0, if so, indicating that the simulation sequence of the jth river structure is determined, setting the return value of Routest (j) as 0, entering step B13, otherwise, indicating that the simulation sequence of the jth river structure cannot be determined, setting the return value of Routest (j) as-1, and entering step B13;

b11, judging whether the diversion unit pointer array of the jth river structure has a value, if so, indicating that the jth river structure is diverted from other river structures, and going to step B12, otherwise, indicating that the jth river structure is not diverted from other river structures, and is the upstream initial river of the current river system, further determining the simulation sequence of the jth river structure, setting the return value of routtest (j) to 0, and going to step B13;

b12, judging whether the simulation sequence of the j river structure shunted by the j river structure is determined in the previous river system traversal process, if so, determining the simulation sequence of the j river structure, setting the return value of RoutetTest (j) to be 0, and entering the step B13; otherwise, the simulation sequence of the jth river structure cannot be determined, the return value of RoutetTest (j) is set to be-1, and the step B13 is carried out;

5. The method for predicting seasonal river interaction with groundwater according to claim 1, wherein the formula for calculating the flow rate between the river and groundwater in step S7 is:

wherein QRIV' represents the flow between the river and the ground water, stored in the field STREAM, whose value being positive represents the river flow to the groundThe leakage amount of the lower aquifer indicates the discharge amount of the underground aquifer flowing to the river when the value is negative; CRIV represents the hydraulic conductivity coefficient between river-underground aquifer, stored in field COND; h is_i,j,kThe calculated water head of a unit where the river reach is located is represented, RBOT represents the elevation of a certain point at the bottom of a river bed, and the elevation is stored in a field STRBTM; HRIV indicates the water level of the river reach and is stored in the field STRSTAGE.

6. The method for predicting seasonal river and groundwater interaction according to claim 5, wherein the HRIV of the water level of the river reach is determined by: