Groundwater Response of Loess Tableland in Northwest China under Irrigation Conditions

: Water induced loess landslides are closely related to the rise of the groundwater level. Therefore, research on the response of the groundwater level to irrigation water holds promise for revealing the mechanism of water-induced loess landslide. Taking Heitai, Gansu Province, as the research area, a coupling model o unsaturated-saturated water movement is established using the HYDRUS-MODFLOW software. The parameters of the model are calibrated and veriﬁed by the Bayesian parameter inversion method combined with ﬁeld observations of the groundwater level. Finally, the change in the groundwater level under di ﬀ erent irrigation amounts is predicted using the optimized model. It is found that a reasonable reduction of the irrigation amount can e ﬀ ectively slow the rise of the groundwater level. This research provides a scientiﬁc reference for the development of reasonable irrigation measures.


Introduction
Loess is a special deposit that is formed in the Quaternary with weak cementation and is widely distributed throughout the world. In recent years, loess landslides triggered by agricultural irrigation have attracted the attention of engineering geologists [1][2][3]. The rise of the groundwater level caused by agricultural irrigation is the essential cause of the development of landslides in the irrigation area [4][5][6]. Therefore, research on the response of groundwater systems under irrigation conditions is the basis for revealing the mechanism of landslides.
Most traditional groundwater models split the complex transformation relationship between soil water and groundwater, resulting in a systematic deviation of the simulation results. To better simulate the flow and storage process of water in a basin, soil water and groundwater should be considered simultaneously to build a comprehensive hydrological model [7][8][9][10][11][12]. At present, many unsaturated-saturated coupling models have been developed, but they are quite different in terms of grid division, coupling mode, and other aspects [13][14][15]. In the coupling of the hydrological model SWAT and the groundwater model MODFLOW, the flow flux between the layers of the unsaturated zone is calculated by the water balance equation [16,17]. SWAT comprehensively considers crop growth, agricultural water management measures, farmland irrigation, and other aspects, but the interaction between the unsaturated system and groundwater system is poor [18]. The MODFLOW-SVAT coupling model couples the soil water movement model SVAT and the groundwater movement model MODFLOW [19]. The unsaturated zone flow movement is calculated by the water balance model. M. Kuznetsov et al. proposed a method to solve three-dimensional variable saturation flow by using the quasi three-dimensional Richards equation and the finite difference scheme. They established the saturated-unsaturated quasi-3D model on this basis [20]. The UZF1-MODFLOW model simplifies the one-dimensional Richards equation as a motion wave to describe the movement process of 1.
Taking Heitai, Gansu Province, as the research area, an unsaturated-saturated coupling flow model was established using the HYDRUS-MODFLOW software combined with rainfall, irrigation, and evaporation data; 2.
Combining groundwater level data monitored in the field with the Bayesian-MCMC random parameter inversion method, the optimization model is obtained by parameter calibration and model verification; 3.
Using the optimized model to predict the change of the trend of the groundwater flow field under different irrigation conditions and exploring effective measures to slow the rise of the groundwater level.

Unsaturated Transport Control Equation
The Richards equation can describe the law of soil water movement, which is as follows: where θ is the volume water content (L 3 L −3 ); h is the pressure head (L); K(h)is the unsaturated hydraulic conductivity (LT −1 ); Z is the vertical depth (L); t is the time (T); S(h) is the source and sink term (L 3 L −3 T −1 ). To solve this equation, we choose the Van Genuchten model to describe the relationship between the water content and pressure head and the Mualem model to describe the relationship between the unsaturated hydraulic conductivity and pressure head. The specific form of van Genuchten-Mualem (VGM) model is as follows: where S e = θ−θ r θ s −θ r is the effective saturation; θ s is the saturated water content; θ r is the residual water content; α is the parameter related to the average particle size; and n is the parameter related to the particle size uniformity, m = 1-1/n.

Basic Theory of the HYDRUS-MODFLOW Model
MODFLOW software has powerful functions to solve various problems related to groundwater flow. The theoretical basis of the numerical simulation of groundwater flow is the three-dimensional movement equation of groundwater in porous media. The equation of motion is derived according to the law of conservation of mass. The three-dimensional partial differential equation of groundwater unsteady flow is as follows: where K xx , K yy and K zz are the components of the hydraulic conductivity in the x, y, and z coordinate directions under the assumption that the main axis direction of the permeability coefficient is consistent with the coordinate axis direction (LT −1 ); H is the groundwater level (L); x, y, and z are the spatial coordinates; W is the source sink term (L 3 L −3 T −1 ); S s is the water storage coefficient, which is dimensionless; and t is the time (T). The most remarkable feature of MODFLOW is its modular structure. MODFLOW first introduced the concept of the stress period, which can be divided into a series of stress periods during the whole simulation period. Each stress period can be further divided into several time steps. MODFLOW uses the centre difference in space and the backward difference in time to obtain a discrete equation [36].
In the coupled HYDRUS-MODFLOW model, two independent control equations are applied to simulate the flow in the saturated zone and the unsaturated zone. The HYDRUS model focuses on the simulation of soil water in the unsaturated zone, and the MODFLOW model mainly simulates the movement process of groundwater in the saturated zone. By embedding the subroutine package HYDRUS into the main program MODFLOW of the groundwater model and taking the water surface as the interface, these equations can effectively simulate the water exchange process between saturated and unsaturated zones. Figure 1 shows a schematic diagram of the HYDRUS-MODFLOW coupling model: The effect of the HYDRUS-MODFLOW coupling model depends to a great extent on how the two models interact in time and space. The simulation area in MODFLOW is divided into several zones, and the average groundwater depth of all the units in a zone determines the pressure head at the bottom boundary of the HYDRUS profile. The calculation method is as follows: where HB is the bottom pressure head of the soil profile (L); ZSURF is the surface elevation (L); DEPTH is the buried depth of the groundwater level (L); and Z1 is the coordinate of the bottom of the soil profile (L). The time step of the HYDRUS model is independent and less than that of MODFLOW. The HYDRUS-MODFLOW coupling model simulates the flow in saturated and unsaturated zones through the discrete interaction of space and time and through the iterative cycle of a stress period.
to simulate the flow in the saturated zone and the unsaturated zone. The HYDRUS model focuses on the simulation of soil water in the unsaturated zone, and the MODFLOW model mainly simulates the movement process of groundwater in the saturated zone. By embedding the subroutine package HYDRUS into the main program MODFLOW of the groundwater model and taking the water surface as the interface, these equations can effectively simulate the water exchange process between saturated and unsaturated zones. Figure 1 shows a schematic diagram of the HYDRUS-MODFLOW coupling model:

Bayesian-MCMC Parameter Inversion Method
The Bayesian method is a recently popular random method, and it is based on the Bayesian principle: where P(m) is the prior distribution of model parameters, which is our understanding of the model parameters before obtaining the observation value d; P (m|d) is the posterior distribution of model parameters, which is our updated and more accurate understanding of the model parameters after obtaining the observation value d; and P(d|m) is the likelihood ratio, which is a measure of the proximity between the model output F (m) and the observation value d. Because most groundwater models are nonlinear, we cannot obtain an analytical expression of the posterior distribution of the model parameters.
The Markov chain Monte Carlo (MCMC) method can be used to sample the posterior distribution of the parameters and estimate the posterior distribution of the model parameters. Vrugt and his partners further proposed the DREAM algorithm, which is a multi-chain MCMC algorithm [37].

Topographic Features
Heitai is located on the north bank of the Yellow River at the intersection of the Yellow River and the Huangshui River ( Figure 2). The total area of the tableland is approximately 10.8 km 2 . The tableland surface is flat and broad, and the slope of the tableland surface is less than 5%. Because of the influence of landslide and collapse, there are many mounds under the tableland, and the ground is uneven (Figure 3).

Stratigraphic Structure and Aquifer Characteristics
The stratigraphic structure of the study area from top to bottom is: Loess, silty clay, pebble layer, and siltstone. According to the characteristics of the lithologic combination of the tableland in the area, the loose and porous loess easily seeps in the vertical direction, the dense silty clay layer is relatively water-proof, the gravel layer has excellent infiltration runoff conditions, and the sand mudstone is relatively water-proof. Therefore, the Heitai tableland is composed of three water-bearing rock groups: loess, sand gravel, and bedrock, from top to bottom. The groundwater in the loess layer is distributed in the Malan loess of the Heitai tableland. The thickness of the aquifer is approximately 20-25 m, which is slightly thick in the middle and east and slightly thin in the west. The buried depth of the groundwater level is 15-26 m, which is generally shallow in the east, shallow in the west, and deep on both sides of the middle. The loess itself has the characteristics of high porosity and developed vertical joints, which makes it a good space for groundwater. The vertical permeability of the loess is far greater than the horizontal permeability, showing the characteristics of homogeneity and anisotropy.

Stratigraphic Structure and Aquifer Characteristics
The stratigraphic structure of the study area from top to bottom is: Loess, silty clay, pebble layer, and siltstone. According to the characteristics of the lithologic combination of the tableland in the area, the loose and porous loess easily seeps in the vertical direction, the dense silty clay layer is relatively water-proof, the gravel layer has excellent infiltration runoff conditions, and the sand mudstone is relatively water-proof. Therefore, the Heitai tableland is composed of three water-bearing rock groups: loess, sand gravel, and bedrock, from top to bottom. The groundwater in the loess layer is distributed in the Malan loess of the Heitai tableland. The thickness of the aquifer is approximately 20-25 m, which is slightly thick in the middle and east and slightly thin in the west. The buried depth of the groundwater level is 15-26 m, which is generally shallow in the east, shallow in the west, and deep on both sides of the middle. The loess itself has the characteristics of high porosity and developed vertical joints, which makes it a good space for groundwater. The vertical permeability of the loess is far greater than the horizontal permeability, showing the characteristics of homogeneity and anisotropy.

Hydrometeorology and Agricultural Irrigation
The study area is located in the northwest inland area, with a semi-arid climate in the middle temperate zone. The supply of precipitation to the groundwater in the area is very limited. In recent years, Heitai has become an important vegetable and fruit planting base. Economic crops that require more water have been mainly planted, and agricultural irrigation has been greatly increased ( Figure 4a). However, irrigation infiltration increases the groundwater supply and destroys the groundwater balance field artificially. Along with the rise of the groundwater level, the thickness of the loess saturated zone increases year by year, which not only causes the overall collapse and landslides on the edge of the tableland but also causes cracks and sinkholes on the surface of the tableland (Figure 4b).
Water 2020, 12, x FOR PEER REVIEW 6 of 15

Hydrometeorology and Agricultural Irrigation
The study area is located in the northwest inland area, with a semi-arid climate in the middle temperate zone. The supply of precipitation to the groundwater in the area is very limited. In recent years, Heitai has become an important vegetable and fruit planting base. Economic crops that require more water have been mainly planted, and agricultural irrigation has been greatly increased ( Figure 4a). However, irrigation infiltration increases the groundwater supply and destroys the groundwater balance field artificially. Along with the rise of the groundwater level, the thickness of the loess saturated zone increases year by year, which not only causes the overall collapse and landslides on the edge of the tableland but also causes cracks and sinkholes on the surface of the tableland (Figure 4b).  The ground elevation, bottom elevation, and initial groundwater level elevation of the

Generalization of the Boundary Conditions of the Calculation Model
The ground elevation, bottom elevation, and initial groundwater level elevation of the simulation area are shown in Figure 5a-c, respectively. There are 5 groundwater level holes (N4, N6, N8, N9, and N10) for monitoring, and their distribution is shown in Figure 5d. The main surface boundary conditions of the study area are rainfall, irrigation, and evaporation. Table 1 shows the statistical data of each influencing factor in a hydrological year, and the irrigation water mainly comes from the Yellow River. The study area is the loess tableland, which is higher than its surroundings and is cut off by valleys. The loess aquifer can neither be supplied by the surface water nor by the lateral runoff of the groundwater outside the area. Some of the phreatic water in the loess layer is generally discharged in the form of a spring on the contact surface of the loess and silty clay, and the other part infiltrates to the sand gravel layer through the relatively impermeable silty clay layer. Therefore, the lateral direction of the model is generalized as the discharge boundary, which is achieved by the Drain module in MODFLOW. Because the permeability of the clay layer is very weak, the bottom of the model is generalized as a water barrier boundary.   For the groundwater flow movement, each month is divided into three stress periods according to the actual situation, and the whole simulation period has 36 stress periods ( Table 2). Each stress period contains several time steps. In the simulation of the HYDRUS module, the initial time step is

Time Division of the Model
A total of 365 days of groundwater dynamic observation data in the simulation area, from 1 July 2016 to 31 June 2017, are selected for model parameter identification. The model is verified using another 365 days of groundwater dynamic observation data from 1 July 2017 to 30 June 2018. For the groundwater flow movement, each month is divided into three stress periods according to the actual situation, and the whole simulation period has 36 stress periods (Table 2). Each stress period contains several time steps. In the simulation of the HYDRUS module, the initial time step is set to 0.01 day, and the minimum and maximum time steps are set to 1 ×10 −3 day and 0.1 day, respectively.

Spatial Division of the Model
In the horizontal direction, the model is divided into 47 rows and 73 columns. A total of 2275 of the 3431 units are active units, and the areas outside the boundary are treated as invalid units. Moreover, as shown in Figure 6, the simulation area is divided into 30 zones according to the surface elevation ( Figure 5a) and the initial groundwater level (Figure 5c) of each unit. That is, the number of one-dimensional unsaturated profiles is 30. In the vertical direction, the HYDRUS profiles are divided into 84-125 finite elements with nodes at different depths.

Spatial Division of the Model
In the horizontal direction, the model is divided into 47 rows and 73 columns. A total of 2275 of the 3431 units are active units, and the areas outside the boundary are treated as invalid units. Moreover, as shown in Figure 6, the simulation area is divided into 30 zones according to the surface elevation ( Figure 5a) and the initial groundwater level (Figure 5c) of each unit. That is, the number of one-dimensional unsaturated profiles is 30. In the vertical direction, the HYDRUS profiles are divided into 84-125 finite elements with nodes at different depths.  Figure 6. Spatial grid division and the HYDRUS module zones. Figure 6. Spatial grid division and the HYDRUS module zones.  Figure 7 shows the changes in the groundwater level at each groundwater level observation point from 1 July 2015 to 11 June 2018. The groundwater level in the Heitai area rises unevenly, and the rising rate in the middle of the tableland is faster than that at the edge of the tableland. Among them, the N6 point has the largest fluctuation amplitude of the groundwater level, which is mainly caused by the infiltration of irrigation water in the surrounding irrigation area. The groundwater level of the N9 also fluctuates obviously. In addition to the impact of vertical irrigation water, N9 is also recharged by the groundwater horizontally due to the terrain of the study area. At the N4 point, the groundwater level fluctuates slightly. The N8 and N10 points have no obvious fluctuation, which shows that they are mainly affected by the lateral recharge of groundwater and that the impact of irrigation recharge is weak.

Dynamic Change of the Groundwater Level
Water 2020, 12, x FOR PEER REVIEW 9 of 15 Figure 7 shows the changes in the groundwater level at each groundwater level observation point from 1 July 2015 to 11 June 2018. The groundwater level in the Heitai area rises unevenly, and the rising rate in the middle of the tableland is faster than that at the edge of the tableland. Among them, the N6 point has the largest fluctuation amplitude of the groundwater level, which is mainly caused by the infiltration of irrigation water in the surrounding irrigation area. The groundwater level of the N9 also fluctuates obviously. In addition to the impact of vertical irrigation water, N9 is also recharged by the groundwater horizontally due to the terrain of the study area. At the N4 point, the groundwater level fluctuates slightly. The N8 and N10 points have no obvious fluctuation, which shows that they are mainly affected by the lateral recharge of groundwater and that the impact of irrigation recharge is weak.

Inversion Results of the Model Parameters
Due to the existence of vertical joints and fissures, loess is generalized as an anisotropic medium in a three-dimensional space, but it is generalized as an isotropic medium in the horizontal direction. The parameters required in the model and their optimized values obtained by the Bayesian-MCMC parameter inversion method are indicated in Table 3. The parameter values after inversion are basically consistent with the actual hydrogeological conditions in the area, and the model has high reliability.

Groundwater Level Simulation Results
Model identification and validation is the key step to evaluate whether the model can truly reflect the hydrogeological conditions of the simulation area. The identification of the model allows repeated adjustment and modification of the hydrogeological parameters, boundary conditions, and source and sink terms of the model so that the model can more accurately reflect the actual hydrogeological conditions in the simulation area. Through the comparison of the actual measured groundwater level of the existing observation holes and the simulated ones, the variation trend and the difference between the two values are minimized. The fitting degree between the measured groundwater data of five observation holes in the study area and the simulated values is taken as the basis of model identification.
As we can see from Figure 8, the groundwater level during the simulation is simulated and compared with the measured values. We take 1 July 2016 to 30 June 2017, a total of 365 days, as the model validation period. The values at the end of every stress period are taken into consideration, and as a result, there are 18 groundwater level observations for each observation point and 90 groundwater level observations for five observation points in total. Figure 8 signifies that the groundwater level simulated by the HYDRUS-MODFLOW model fits well with the measured value. The evaluation index of the simulation effect RMSE = 0.34 m, so the simulation result of HYDRUS-MODFLOW at the regional scale is reliable, which indicates that the optimized model reflects the real dynamic change characteristics of the groundwater level. The observed values of the groundwater level are basically included in the 95% confidence interval, which indicates that the simulation uncertainty ranges are statistically adequate [37].

Inversion Results of the Model Parameters
Due to the existence of vertical joints and fissures, loess is generalized as an anisotropic medium in a three-dimensional space, but it is generalized as an isotropic medium in the horizontal direction. The parameters required in the model and their optimized values obtained by the Bayesian-MCMC parameter inversion method are indicated in Table 3. The parameter values after inversion are basically consistent with the actual hydrogeological conditions in the area, and the model has high reliability.

Groundwater Level Simulation Results
Model identification and validation is the key step to evaluate whether the model can truly reflect the hydrogeological conditions of the simulation area. The identification of the model allows repeated adjustment and modification of the hydrogeological parameters, boundary conditions, and source and sink terms of the model so that the model can more accurately reflect the actual hydrogeological conditions in the simulation area. Through the comparison of the actual measured groundwater level of the existing observation holes and the simulated ones, the variation trend and the difference between the two values are minimized. The fitting degree between the measured groundwater data of five observation holes in the study area and the simulated values is taken as the basis of model identification.
As we can see from Figure 8, the groundwater level during the simulation is simulated and compared with the measured values. We take 1 July 2016 to 30 June 2017, a total of 365 days, as the model validation period. The values at the end of every stress period are taken into consideration, and as a result, there are 18 groundwater level observations for each observation point and 90 groundwater level observations for five observation points in total. Figure 8 signifies that the groundwater level simulated by the HYDRUS-MODFLOW model fits well with the measured value. The evaluation index of the simulation effect RMSE = 0.34 m, so the simulation result of HYDRUS-MODFLOW at the regional scale is reliable, which indicates that the optimized model reflects the real dynamic change characteristics of the groundwater level. The observed values of the groundwater level are basically included in the 95% confidence interval, which indicates that the simulation uncertainty ranges are statistically adequate [37].

Model Validation
Before using the established model for simulation and prediction, we should first verify the model. We take 1 July 2017 to 30 June 2018 as the model validation period, a total of 365 days, to further verify the parameters of the optimized model. It can be seen from Figure 9 that the fitting effect between the calculated groundwater at each observation hole and the measured values during the model validation period is good. RMSE = 0.33 m indicates that the trend of the numerical simulation values of the Loess phreatic water level is close enough to the measured values. The model can accurately reflect the hydrogeological and dynamic characteristics of the simulation area from the perspective of the overall simulation situation and can be used for water level prediction. The simulation results show that the HYDRUS-MODFLOW model has a good simulation effect on the water exchange process between saturated and unsaturated zones at the regional scale, which is consistent with Reference [24].
Water 2020, 12, x FOR PEER REVIEW 11 of 15

Model Validation
Before using the established model for simulation and prediction, we should first verify the model. We take 1 July 2017 to 30 June 2018 as the model validation period, a total of 365 days, to further verify the parameters of the optimized model. It can be seen from Figure 9 that the fitting effect between the calculated groundwater at each observation hole and the measured values during the model validation period is good. RMSE = 0.33 m indicates that the trend of the numerical simulation values of the Loess phreatic water level is close enough to the measured values. The model can accurately reflect the hydrogeological and dynamic characteristics of the simulation area from the perspective of the overall simulation situation and can be used for water level prediction. The simulation results show that the HYDRUS-MODFLOW model has a good simulation effect on the water exchange process between saturated and unsaturated zones at the regional scale, which is consistent with Reference [24].

Prediction of the Groundwater Level Change Trend
Using the above optimized model, we can predict the change of the groundwater level under different irrigation conditions. In this paper, the change of the groundwater level in the study area over 3 years is predicted under the condition that the irrigation amount is 1 times and 0.8 times the initial irrigation amount. Figure 10 shows the distribution of the groundwater flow field in the first year and the third year under the conditions of different irrigation volumes.

Prediction of the Groundwater Level Change Trend
Using the above optimized model, we can predict the change of the groundwater level under different irrigation conditions. In this paper, the change of the groundwater level in the study area over 3 years is predicted under the condition that the irrigation amount is 1 times and 0.8 times the initial irrigation amount. Figure 10 shows the distribution of the groundwater flow field in the first year and the third year under the conditions of different irrigation volumes.
It can be seen from Figure 10 that under the condition of the original irrigation amount and 0.8 times the original irrigation amount, the contour of the groundwater level obviously shifts right, that is, the groundwater level rises with time. Among these, the contour of the 1702 m groundwater level on the left side expands with the increase of the simulation time. However, under the condition of the original irrigation, the movement of the groundwater level contour is more obvious. The results fully show that a reasonable reduction of irrigation can effectively slow down the continuous rise of the groundwater level.
Under the influence of long-term irrigation, the groundwater level of Heitai is rising year by year, which leads to many landslides along the border of the tableland [4]. Treatments such as slope cutting and load reduction have been carried out in this area and potential deformation bodies have been temporarily removed. However, because the hydrogeological conditions have not fundamentally changed, the slope body is likely to undergo a potential deformation at any time. Through the discussion in the previous section, it is concluded that a reasonable reduction of irrigation can effectively slow down the continuous rise of the groundwater level. There are many ways to reduce the amount of irrigation: first, we should change the method of irrigation, reducing flood irrigation to prevent water waste; second, we should improve the type of crops planted, planting more types of crops with strong drought resistance and less water demand. Simultaneously while reducing agricultural irrigation, other measures should also be taken to help reduce the groundwater level. Direct groundwater drainage schemes are the most effective and direct measures for landslide disaster control. It can be seen from Figure 10 that under the condition of the original irrigation amount and 0.8 times the original irrigation amount, the contour of the groundwater level obviously shifts right, that is, the groundwater level rises with time. Among these, the contour of the 1702 m groundwater level on the left side expands with the increase of the simulation time. However, under the condition of the original irrigation, the movement of the groundwater level contour is more obvious. The results fully show that a reasonable reduction of irrigation can effectively slow down the continuous rise of the groundwater level.
Under the influence of long-term irrigation, the groundwater level of Heitai is rising year by year, which leads to many landslides along the border of the tableland [4]. Treatments such as slope cutting and load reduction have been carried out in this area and potential deformation bodies have been temporarily removed. However, because the hydrogeological conditions have not fundamentally changed, the slope body is likely to undergo a potential deformation at any time. Through the discussion in the previous section, it is concluded that a reasonable reduction of irrigation can effectively slow down the continuous rise of the groundwater level. There are many ways to reduce the amount of irrigation: first, we should change the method of irrigation, reducing flood irrigation to prevent water waste; second, we should improve the type of crops planted, planting more types of crops with strong drought resistance and less water demand. Simultaneously while reducing agricultural irrigation, other measures should also be taken to help reduce the groundwater level. Direct groundwater drainage schemes are the most effective and direct measures for landslide disaster control.

Conclusions
The main conclusions of this study are:

1.
On the basis of the coupling principle and operation mechanism of the HYDRUS-MODFLOW coupling model, the model is applied to simulate the groundwater level of Heitai in Gansu Province. The simulation results show that the HYDRUS-MODFLOW model has a good simulation effect on the water exchange process between saturated and unsaturated zones at the regional scale.

2.
To further improve the practicability and simulation accuracy of the HYDRUS-MODFLOW model, it is combined with the Bayesian-MCMC parameter inversion method. The parameters in the model are inverted and verified using the measured groundwater level data in the field water level holes. The results show that the simulation values of the coupling model fit well with the measured values, which indicates that the model can better simulate the transformation relationship among surface water, soil water, and groundwater at the regional scale in Heitai.

3.
The development trend of the groundwater level of the Heitai groundwater system in Gansu Province in the next 3 years under different irrigation intensities is predicted using the optimized model. The prediction results show that the groundwater level is seriously affected by the irrigation intensity. The groundwater level increases with the increase of the irrigation intensity and decreases with the decrease of the irrigation intensity. A reasonable reduction of the irrigation intensity can slow the rising speed of the groundwater level.

4.
Measures to reduce the groundwater level in Heitai are recommended, such as reasonable reduction of the irrigation amount by changing the irrigation mode; adjustment of the crop structure and planting area to reduce uneven irrigation as much as possible; and direct discharge of groundwater by adopting the drainage test scheme.