Interaction Simulation of Vadose Zone Water and Groundwater in Cele Oasis: Assessment of the Impact of Agricultural Intensification, Northwestern China

: Agricultural intensification has boosted land productivity, but it has also created new sustainability issues. As one of the most important human habitations and agricultural farming areas in arid areas, the Cele Oasis has a very developed agricultural system. This paper studies the long-term effects of different types of agricultural intensification strategies on groundwater level fluctuations in the Cele Oasis. A soil water flow (HYDRUS-1D) and aquifer simulation (MODFLOW) coupling model were used to construct the geometric structures of the vadose zone and saturated zone in the Cele Oasis and to analyze the recharge and discharge mechanism of the oasis. The results showed that HYDRUS-1D accurately simulated soil moisture transport in the Cele Oasis, providing reliable data for calibration of the MODFLOW model. The groundwater level simulated by MODFLOW was in good agreement with the observed value. The results of the R 2 , RMSE, and NSE were ranges of 0.77–0.90, 0.45–0.74 m, and 0.76–0.87, respectively. The errors were acceptable limits. The coupling model predicted the responses of different agricultural types and cropping scenarios to groundwater. Predictions showed that the groundwater level in the Cele Oasis remained stable under the current cropping scenario (100% cropping intensity), and that the groundwater level decreased slightly under the cropping scenario (110% cropping intensity and 130% cropping inten-sity). When the cropping scenario was at 170% cropping intensity, the groundwater level decreased rapidly, and the maximum drawdown value was 7 m. Therefore, the maximum cropping intensity of the Cele Oasis in the future should be 130 % of the current cropping intensity.


Introduction
In arid and semi-arid regions, where surface water is greatly affected by short-term changes in weather patterns, sustainably obtained groundwater is critical for sustaining human life, enhancing the intensification and diversification of agricultural systems, combating climate change and variability, and promoting economic development [1]. Groundwater is often considered a stable and reliable source of fresh water. Globally, 301 million hectares of irrigated land exist, 38% of which is irrigated by groundwater [2,3]. Global groundwater extraction has more than quadrupled over the past 50 years as a result of population growth, industrial centralization, agricultural intensification, and increasing numbers of efficient water pumps. Groundwater systems are under increasing pressure in different regions, with more water being withdrawn from aquifers than is replenished. This has resulted in a gradual decline in groundwater levels over time [4,5]. In addition, future climate change and population growth will force governments to make significant efforts to attain food security, and the impact of intensified agriculture on groundwater systems will, therefore, become more pronounced in drylands [1,6,7].
The sustainable development of aquifers is generally understood as a dynamic equilibrium between extraction and recharge, which means that groundwater is considered a renewable resource in nature [8]. However, this concept oversimplifies the hydrogeological system and the vadose zone between the surface and the groundwater table affects groundwater transport, which determines the availability of groundwater recharge [9,10]. Recharge estimation plays an important role in the management of groundwater systems, and the amount and timing of recharge are affected by climatic and geological factors [11][12][13]. Therefore, modeling vadose water transport in aquifers allows for more realistic assessments of groundwater recharge and accurate predictions of groundwater level rise or fall [14]. Vadose zone flow simulation is a complex process that requires the necessary data. Numerical simulation is an effective method for evaluating groundwater recharge [15]. HYDRUS-1D [16,17], which is an advanced hydrological model, can simulate a onedimensional variable-saturated soil water, heat, and solute transport. HYDRUS-1D is also widely used to simulate water uptake by crop roots in soil aquifers [18][19][20][21][22].
The hydrological processes that take place in the vadose zone and saturation zone are closely related. The error caused by simplifying estimates and the uncertainty of recharge in vadose zones is directly transferred to the predicted groundwater regression. If this error is ignored, it leads to deviations in the estimated return of irrigated groundwater [23,24]. Therefore, recharge flux simulated using the independent vadose zone model HY-DRUS-1D can significantly improve the accuracy of the saturation zone model in predicting groundwater level [1]. MODFLOW and its variants are one of the most popular types of groundwater modeling packages for solving both confined and unconfined equations [25]; a major disadvantage of MODFLOW, however, is that it simplifies the calculation of water transport in the vadose zone. Many attempts have been made to incorporate a vadose zone flow model into MODFLOW and have achieved good results [26][27][28]. The coupled model effectively simulates the vadose-unsaturated hydrological process on a regional scale [14] and can be used to predict the long-term impact of changes in surface water management on groundwater level.
The desert oasis is a unique intrazonal landscape, supporting the survival and development of vegetation and human beings in arid or semi-arid areas [29,30]. Oases account for 4%-5% of arid areas in China, but support more than 90% of the population in arid areas and create more than 95% of the social wealth [31]. Although glacier melt water is the main water source for agriculture in summer (May-September), groundwater is the main water source in autumn (October-November) and winter (December-April). The increase in groundwater pumping in autumn and winter leads to a serious decline in groundwater storage in an oasis [32]. In the Cele Oasis, adjustments to the cropping structure and the promotion of a "fruit trees-crops intercropping" cropping structure improve the irrigation system [33]. The irrigation water demand of crops in autumn and winter on groundwater resources directly affects the sustainable development of oasis groundwater systems.
Irrigated agriculture is not only an important part of oasis agricultural economic development, but also one of the main drivers of groundwater resource reduction [34]. Previous research on oases in arid regions of the Northwest China [35][36][37][38] has used mainly MODFLOW to simulate groundwater movement, while ignoring vadose zone water transport. The government has increased cropping intensity in order to achieve economic growth and to meet the challenges of food security in the future. This will pose a significant threat to groundwater security in the Cele Oasis, thus requiring scientific research to mitigate it. The purpose of this study is to evaluate the impact of agricultural intensification on the groundwater level of the Cele Oasis. We used HYDRUS-1D and MODFLOW models to simulate the effects of different cropping intensities on the water table. The results will help in the development of robust aquifer extraction schemes to address the competition between agricultural intensification and groundwater security.

Study Area
The Cele Oasis lies at coordinates 36°54′-37°09′ N, 80°37′-80°59′ E: the total area between the main part of the oasis and the desert transition zone is 189.64 km 2 [30], and the cultivated land area is 22,974.75 ha. The oasis is located on the southern edge of the Taklimakan Desert and the northern foot of the Kunlun Mountains, China; it has a typical continental arid desert climate and an average altitude of 3200 m in the mountains and 1500-1800 m in the plains.
The temperature here is largely different between day and night. The average annual precipitation is 35.5 mm, and the annual evaporation is 2751 mm. The annual average evaporation data come from the 20 cm evaporation pan data of the Cele Station from 1960 to 2010 [39]. The vegetation coverage is low and the ecological environment is very fragile. The Cele River provides a source of surface water for the evolution of the oasis, with an average annual runoff of 1.27 × 10 8 m 3 , which flows into the extremely arid Taklimakan Desert.
The land use in the Cele Oasis is predominantly cropland, forest, and build. The cropland landscape covers about 9.6% of the total study area, forest covers 87.8%, and build land covers 13%. The main types of crops in Cele Oasis are red dates, walnuts, wheat, and corn. Crops typically are irrigated five or six times during the growth period, mostly from groundwater. (Figure 1).
The canal system inside the oasis is developed, forming a typical artificial farmland oasis landscape. The direction of flow of the groundwater aquifer in the oasis is the same as that of the surface flow. The buried depth of groundwater is less than 50 m, and the aquifer is a quaternary pore phreatic water system. The soil type in the study area is loamy sand, which has a high infiltration rate. In recent years, due to the impact of climate change and human activity, the agricultural cropping structure of the oasis has changed greatly, which has greatly affected the safety of the aquifer.

Data Sources
In this study, meteorological data from 2005-2017 were collected from the Cele National Station, taking into account the time-matching of groundwater data and meteorological data and using the FAO-56 Penman-Monteith method [40] to calculate the daily reference evapotranspiration. GF-1 remote sensing data (with a resolution of 2 m) were obtained from the China Resources Satellite Research Center (http://www.cresda.com/CN/ (accessed on 14 July 2016)). ENVI5.2 and ArcGIS software were used to preprocess the data. According to the land use/land cover classification system of the Chinese Academy of Sciences, a land use type map of the Cele Oasis was drawn using the supervised classification method. Fifty sample sites were randomly selected to test the accuracy of the classification results: the kappa coefficient was 0.92, which meets the standards of research accuracy. The Harmonized World Soil Datasets, version 1. After analyzing the pumping test data, the hydraulic conductivity and specific yield of the aquifer in the region were obtained. From 2015 to 2016, the soil moisture was monitored using soil temperature hygrometry at the field test site of the Cele National Station. The equipment provided automatic monitoring and the monitoring interval was set to 1 h. A hydrogeological map of the study area was obtained from the groundwater resource atlas of Xinjiang.

Modelling Description
HYDRUS-1D software [41] simulated the one-dimensional movement of water in variable-saturated porous media, soil moisture storage, and root water uptake. The model used a finite element numerical method to solve Richard's equation of variable-saturated flow. In order to input the evapotranspiration (ET) and groundwater recharge calculated by HYDRUS-1D into the MODFLOW (cell) for data exchange, the study area was divided into different sites called "land units" according to the land use map and soil type map, and the recharge of different land type units was simulated. The main land use types in the oasis were cropland (9.6%), forest (87.7%), and build (2.7%). Loamy sand was the main soil type in the Cele Oasis. According to the land use type and soil type spatial distribution map, 3 simulations were constructed using HYDRUS-1D software for cropland loamy sand (CLS), forest loamy sand (FLS), and build loamy sand (BLS) ( Table 1). Groundwater modeling system software (GMS), a type of modular three-dimensional steady-state finite difference software, was used to quantify the groundwater in the Cele Oasis. The coupling modeling method of HYDRUS-1D and MODFLOW software used the drainage in the HYDRUS-1D model as the recharge in the MODFLOW model. The bottom flux calculated using the HYDRUS-1D model subtracted the pumping water from the well to obtain the "net recharge" of groundwater. The bottom boundary condition of the pressure head designated as HYDRUS-1D was used as the groundwater level in MODFLOW [42] to evaluate the influence of the vadose zone on groundwater fluctuation. Due to the deep groundwater in the arid area, the groundwater recharge of soil water in the unsaturated zone can be ignored for groundwater depths greater than 5 m. Therefore, we only considered the one-way process of groundwater recharge from the vadose zone in the Cele Oasis. The MODFLOW model divided the aquifer into blocks or grids and divided the soil or hydrogeological parameters into a region with similar parameters. Different regions corresponded to different HYDRUS-1D soil profiles. Figure 2 shows the coupling model of water flow through soil profile migration. Figure 3 shows a flow chart of the aquifer simulation using HYDRUS-1D and MODFLOW.

HYDRUS-1D Model
The HYDRUS-1D model is a one-dimensional vertical model, which simulates water, heat, solute transport, and root water uptake in vadose zones [43]. The model used Richard's equation with added source-sink terms to solve for vertical soil moisture transport [44] (Equation (1)). The model assumed that air does not hinder the flow of water; in addition, the model used the van Genuchten equation [45] or the Brooks and Corey characteristic equation [46] to calculate the unsaturated soil hydraulic characteristics. Table 2 gives the soil hydraulic parameters of three land use types in the Cele Oasis. The onedimensional Richard's equation was given by: where is the volumetric water content (dimensionless), K is the hydraulic conductivity The water content, (h), and the unsaturated hydraulic conductivity, K(h), made Richard's equation a highly nonlinear equation that needed to be solved numerically. Note: θr is the soil residual water content, θs is the soil saturated water content, α and n are curve fitting parameters, Ks is the saturated hydraulic conductivity, and l is the molecular diffusion coefficient.

Evaporation and Transpiration
According to climate and soil parameters, ET0 was calculated using the Penman-Monteith model [47,48]. The potential evapotranspiration of crops was calculated using Eq (2): where is the crop coefficient. According to FAO56, the values of in this paper were 0.4, 0.8, and 0.5, in the early stage (ground coverage is less than 10%), the middle stage (ground coverage is from 10% to maturity), and in the later stage (ground coverage is from maturity to harvest), respectively [40].
was the reference evapotranspiration. and precipitation from 2005 to 2017 are shown in Figure 4.
The HYDRUS-1D model calculated the soil evapotranspiration (Ep) and crop evapotranspiration (Tp) separately, given in Equations 3 and 4, respectively: where LAI is the leaf area index, and k is the canopy extinction coefficient.

Soil Hydraulic Parameters
From 2015 to 2016, a 2-m soil moisture monitoring experiment was carried out at the Cele National Station. Therefore, HYDRUS-1D used 2 m as the soil profile depth to simulate the bottom flux of each land unit and considered the soil texture in the Cele Oasis to be uniform; the hydraulic characteristics of the whole soil profile were also considered to be uniform. In the HYDRUS-1D model, the soil volume mass and soil particle size classification were the inputs, and parameters such as saturated water content θs, residual water content θr, and saturated hydraulic conductivity Ks were preliminarily determined using the Rosetta module of a neural network. Compared with the observed water data, the error between the simulated value and the measured value was reduced by adjusting the parameters. The results are shown in Table 2.

Root Water Uptake
HYDRUS-1D used a water stress response function to study root water uptake, and the Feddes model was used to calculate root water uptake [49]. The parameters in the root water uptake-water stress response function were taken from the HYDRUS-1D database [50].

Initial and Boundary Conditions
The initial condition of the HYDRUS-1D model was the initial measured value of soil moisture. The upper boundary used the atmospheric boundary of the surface integrated water, and meteorological data were used to specify the atmospheric boundary conditions. The lower boundary was located in the unsaturated zone because it did not reach the phreatic layer (the groundwater table depth is greater than 5 m), and it was set as the free drainage boundary.

MODFLOW Model
MODFLOW is one of the most widely used groundwater numerical simulation models and is a physical model. MODFLOW was combined with Darcy's law and the groundwater flow equilibrium equation, using a finite difference method, to solve the groundwater motion equation. Using the MODFLOW module in GMS software, this paper constructed a MODFLOW model of the Cele Oasis [51] and solved the three-dimensional groundwater flow equation (Equation (5)) to obtain the head value.
where , , and are hydraulic conductivity coefficients (L/T) in the x-, y-, and zdirections, respectively; h is the hydraulic head (L); is the specific storage (1/L); and W is the recharge/discharge rate per unit volume (1/T).
In the MODFLOW finite difference method, continuous time and space were divided into a series of discrete blocks or grids, including horizontal and vertical subdivision grids. On the horizontal plane, the Cele Oasis was composed of 24,585 elements on a 200 m × 200 m grid, including 149 rows and 156 columns. The vertical profile of the aquifer was that of a single phreatic aquifer, and the simulation depth was (5-100 m), according to the groundwater level monitoring data and evaporation, and other related hydrological data. From 1 May 2005 to 30 April 2017, the time length was 4382 days, which was divided into 24 pressure periods. The pressure period was divided into a flood stress period (May-September) and a dry season stress period (October-April). The time step was one day. Therefore, the flood stress period was 153 days, and the dry season stress period was 212 days.

Initial and Boundary Conditions
According to the groundwater monitoring data of Cele Oasis in May 2005, the kriging interpolation method was used to interpolate 25 observations of well data in ArcGIS software to obtain the initial head. The upper boundary of the oasis is the phreatic surface, where the water exchange boundary with changing positions occurs, including irrigation infiltration and pumping wells. The lower boundary is bounded by a weak permeable layer, which is generalized as an impermeable boundary. The southern mountain outlet area is the flow boundary of lateral inflow. The eastern, northwestern, and northern areas are bound by underground watersheds, which are generalized as zero-flow boundaries. A small amount of lateral outflow is observed in some areas for lateral outflow boundary treatment.

Model Calibration
Observed data collected for 423 days, from May 4, 2016 to April 30, 2017, were selected to verify the HYDRUS-1D model. During this period, a total of 10 irrigation times were used. The observed values were compared with simulated values to verify the accuracy of the HYDRUS-1D model simulation. In the case of MODFLOW, the groundwater monitoring data of 25 observation wells, collected from 2005 to 2017, were used for model calibration, and the simulated and observed groundwater levels as functions of time were compared. The commonly used root mean square error (RMSE) [52] and the Nash-Sutcliffe model efficiency coefficient (NSE) [53] were used to evaluate the accuracy of the simulation results of the model. The RMSE and NSE are given by Equations (6) and (7), respectively.
where and are the ith values of the observed and the simulated dataset, respectively, Ō is the average of the observed values, and n is the total number of observations.

Scenario Simulation
The industrial foundation of the Cele Oasis is weak, and the low-value-added crops grown there are mainly walnuts, red dates, and wheat. In order to promote the development of the Cele Oasis, the state implemented a policy of "targeted poverty alleviation." Many institutions promote the improvement of cropping intensity, such as through "fruit trees-crops intercropping", to drive agricultural production and improve residents' income. However, this kind of "fruit trees-crops intercropping" increases the groundwater pumping capacity in the dry season (October-April). The long-term impact of this highintensity groundwater extraction strategy during the drought period is still unclear. Therefore, this paper used 4 different cropping scenarios to evaluate the impact on groundwater. The first scenario was based on the current 100% cropping intensity. In the second, third, and fourth cropping scenarios, 105CI, 110CI, and 115CI were considered, respectively. The four different cropping scenarios assumed that the planting area of forest fruit crops remained unchanged: only the cropping area of wheat/maize changed. (Table 3). The CropWat Model established by the Land and Water Development Division of FAO was used to calculate the increased irrigation water demand due to the expansion of planting areas. According to the irrigation water storage for different crops obtained from the experimental plot of the Cele Station, the monthly average irrigation water demand for forestry fruit crops was 1.5 × 10 3 m 3 /hm 2 , and the monthly average irrigation water demand of wheat and maize was 1.05 × 10 m 3 /hm 2 . Equation 8 is used to calculate the water extraction of crops in the dry season for different cropping intensities: where is the planting area in scenario i, is the total area (ha), ( ) is the sum of the crop irrigation water requirement in the dry season (mm), 242 is the number of drafting water days in the dry season, and is the increased draft (mm/day) in scenario i.

HYDRUS-1D Calibration Results and Recharge Flux
The HYDRUS-1D model was used to simulate the cropland and forest land units in the study area. The spatial location of soil water monitoring sites is shown in Figure 1 Figure 5 shows that the observed water content and the model simulation results fit each other well. The soil water content increased rapidly after irrigation at 50 cm, whereas soil water content at 150 cm responded only two days after irrigation as the rate of change range was relatively small. Due to the complex soil environment, the NSE values between the simulated and observed values for soil water content in the shallow layers were 0.823 and 0.792, respectively, whereas the NSE values for the deep soil water content were 0.876 and 0.81. The fitting of the simulation values with the observed values was good. In terms of the RMSE, the simulation results for deep soil were better than those for shallow soil. In general, the HYDRUS-1D model had high simulation accuracy for soil moisture content. This model could, therefore, readily be used in the analysis of cropland and forest water migration and in the calculation of cropland and forest water supply to groundwater.
The simulation results showed that the largest cumulative recharge flux came from forest land (walnut, red date, pomegranate) (Table 4), mainly from irrigation infiltration recharge, as high evaporation loss in the Cele Oasis and field bunds prevented runoff from cropland and forest land. Compared to that of cropland and forest land units, the cumulative recharge flux from build land units was low. This is mainly due to hardened roads and a lack of precipitation recharge. Spatially, the bottom fluxes of different land units were fed into each model unit in MODFLOW.

MODFLOW Calibration Results
The model correction needed to reflect the real groundwater level value. Therefore, the real groundwater level value of the observation well was input into the model, and the actual discharge of the pumping well was adjusted continuously; thus the simulated water level value was basically consistent with the water level value of 25 observation wells. The adjusted K value is shown in Figure 6a. The comparison between the calculated water level after model correction and the observed water level is shown in Figure 6b. Groundwater monitoring data from May 2005 to April 2017 were used to correct the simulated groundwater recharge in the cropping scenario (100CI). The observed data from 5 observation wells in the study area were compared with MODFLOW simulation data ( Figure 6a). Figure 7 shows the observed and simulated head values of 5 representative observation wells at 1:1. Figure 8 shows the relationship between the simulated and observed average groundwater levels. As shown in Table 5, the ranges of R 2 and RMSE, between the observed and simulated values, were 0.77-0.90 and 0.45-0.74 m, respectively. The lowest NSE value of well-3 and well-4 was 0.76, and the RMSE value of well-1 was highest at 0.74, indicating that the groundwater level simulated by the HYDRUS-1D and MODFLOW coupling models was consistent with the observed values. The simulation accuracy was high. Using GMS software to calculate the water budget of the Cele Oasis aquifer, the results showed that the recharge of the Cele Oasis water resources mainly consisted of melting water from alpine snow and ice, accounting for 96% of the total recharge. During the simulation period, the long-term average annual net inflow volume of the oasis was 1.3 × 10 8 m 3 , whereas the long-term average annual net outflow volume was 1.33 × 10 8 m 3 . This indicated that the Cele Oasis maintained a dynamic water balance under cropping intensity (100CI).

Effects of Different Cropping Scenarios on Groundwater Level
In order to examine the impacts of different cropping scenarios on the groundwater level, the average groundwater level of the flood period from May 2005 to April 2017 was used as the base year; under four different cropping scenarios, the groundwater level predicted by the flood period in April 2017 was compared with the base year to assess the impact of each cropping scenario on the groundwater level of the aquifer.
As compared to the base year, the spatial distribution of the groundwater level under different cropping scenarios is shown in Figure 9. Ten representative sites were selected (Figure 9a) to evaluate the fluctuation of groundwater level under different cropping scenarios (Table 6). Under the current cropping scenario, the groundwater level in 2017 was 0.23, 2.98 m lower than that in the base year. In the area where the water level was greater than 20 m, the surface water recharge was weak, which caused the water level fluctuation. In the area where the water level was less than 20 m, the surface water can effectively recharge the groundwater and make the groundwater level change only a little, which can be ignored. This indicates that there was considerable interaction between surface water and shallow groundwater. In general, the groundwater level in the Cele Oasis did not show a significant upward/downward trend from 2005 to 2017. This indicates that, under the current cropping intensity (100CI), the Cele Oasis maintains a dynamic water balance and that the sustainable development of the oasis remains stable. Figure 9a shows the spatial distribution of groundwater level in the oasis under the 100CI scenario.
The MODFLOW model was used to simulate groundwater level change for different cropping intensities in the Cele Oasis. The results are shown in Figure 9b-d and Table 6. With increasing cropping intensity (110CI, 130CI, and 170CI), the groundwater level in 2017 was significantly lower than that in the base year (2005). When the cropping intensity was 110%, 130%, and 170%, the groundwater level was 1.32-4.05 m, 2-4.54 m, and 3.09-7.07 m lower than that in the base year, respectively. Figure 9 shows that when the cropping intensity was greater than 100%, the areas with a groundwater level less than 5 m disappeared. When the cropping intensity was 110% and 130%, the groundwater level does not have a significant decline. When the cropping intensity was 170%, the groundwater level of the oasis dropped rapidly and the area with a water level of 5-10 m was reduced by 208%, compared with the same water level in the 100CI scenario.

Modeling Method
In this paper, HYDRUS-1D and MODFLOW models were used to evaluate the influence of long-term groundwater extraction from the Cele Oasis. The HYDRUS-1D model was used to simulate the soil water flow in cropland, forest land, and build land units. Due to the existence of the vadose zone, the exchange between the surface water and groundwater decreased [54]. The difficulties in numerical simulation of coupled flow systems in the vadose zone-saturated zone lie in the nonlinearity of the soil water model, the inhomogeneity of parameters, and the variability of hydrological stress [55][56][57]. These difficulties make the recharge of groundwater from soil water challenging to tackle in MOD-FLOW modeling [58], and sometimes time-consuming [59]. Mail [1] and Li et al. [59] demonstrated that HYDRUS-1D can accurately estimate groundwater recharge under different cropping conditions for different study areas. In this study, the recharge calculated by HYDRUS-1D for different land units was accurate and was close to the recharge calculated by Qian et al. [60]. HYDRUS-1D simulated soil water content at different depths and the NSE of the observed value was 0.82. This means that there was good agreement between the simulated value and the observed value.
The Cele Oasis is a typical piedmont alluvial plain. The influence of the deep groundwater level on the soil water in the vadose zone can be ignored, which makes the coupling of the vadose zone and the saturated zone become a single coupling. The bottom flux calculated from the HYDRUS-1D model subtracts the pumping water to obtain the "net recharge" of groundwater. The bottom boundary condition of the pressure head designated as HYDRUS-1D was used as the burial depth of the groundwater level in MOD-FLOW [42] to evaluate the influence of the vadose zone on groundwater fluctuation. The MODFLOW model simulation results were consistent with those of Liu et al. [61], who used an optimization algorithm to calculate groundwater level in the Cele Oasis.

Increased Cropping Intensity: Enlightenment and Strategies
In arid and semi-arid regions, oases play a vital role in human survival and development. The main industry in oases, agriculture, places a very high demand on groundwater compared with agriculture in other regions [62]. Due to the implementation of the "targeted poverty alleviation" policy, how to increase the cropping intensity and economic development in oases with minimal damage to their sustainable development has become paramount. The results also show that the Cele Oasis needs to maintain a long-term balance between water consumption and water availability, while addressing the shortage of water resources caused by increased cropping intensity.
With improvements in drilling technology and the popularization of high-power pumping pumps, unplanned groundwater extraction causes serious waste of groundwater resources. It leads to an unsustainable groundwater system in the short term [63]. Groundwater management, especially groundwater management in oases, is an increasingly important issue. In this context, the results of this study could be used to optimize the layout of pumping wells in oases, so as to realize the sustainable utilization of groundwater in irrigation areas [64,65].
The Cele River is the main source of surface water in the Cele Oasis, with an annual average net inflow volume of 1.2 × 10 8 m 3 . However, more than 59.5% of the irrigation water in the region is from groundwater [60]. The results of this study show that increasing cropping intensity will lead to the unsustainable development of groundwater in the oasis, and that the combined use of surface water and groundwater will help to alleviate shortages of groundwater in the oasis. Liu [61] and Ma et al. [66] improved the layout of pumping wells using optimization algorithms, but there are still pumping wells being overused, leading to serious groundwater extraction. This result predicts that increases in cropping intensity in the future will lead to a significant decline in groundwater levels. Therefore, it is necessary to allocate surface water, improve water conveyance networks, and boost irrigation efficiency so as to ensure the sustainable development of groundwater resources. The average annual runoff of the Cele River is 1.2 × 10 8 m 3 , but a large amount of the surface water flows out of the Cele Oasis each year during the flood period which could be used as a potential source of groundwater recharge. Therefore, it is necessary to construct water conservancy projects in the Cele Oasis for groundwater recharge and agricultural irrigation.

Conclusions
HYDRUS-1D and MODFLOW models were used to evaluate the effects of different cropping intensities on groundwater in the Cele Oasis. HYDRUS-1D was used to simulate soil moisture migration in the vadose zone, and the cumulative bottom fluxes of different land use units (cropland, forest, build) from 2005 to 2017 were estimated. The results showed that the bottom fluxes of forest land were higher than those of cropland, and the bottom fluxes of build land were the lowest. The bottom fluxes of the three types of units were 43.2 mm, 159.6 mm, and 4.8 mm, respectively. The bottom fluxes of the three types of units were used as a recharge source for MODFLOW. By comparing the observed groundwater level data with the simulated data, it was found that MODFLOW accurately simulated groundwater level change. This indicates that the HYDRUS-1D and MOD-FLOW models successfully simulated the interaction between soil water and groundwater in the Cele Oasis. By comparing the effects of different cropping intensities on groundwater, the results showed that there was no significant change in the groundwater level under the current cropping intensity (100CI) scenario. When the cropping intensity increased to 110% and 130%, the groundwater level decreased slightly, and the maximum drawdown depth was 1.56 m. When the cropping intensity reached 170%, the groundwater level decreased significantly. The results indicated that the maximum degree of agricultural intensification in the Cele Oasis was 130% of the current cropping intensity. When it was greater than 130%, the groundwater in the Cele Oasis was not compensated by irrigation water, resulting in a sharp drop in the water level and leading to concerns about the groundwater safety and stability of the oasis. The spatial distribution of the groundwater level in the oasis shows that the northwest region is more vulnerable to increases in groundwater exploitation. Therefore, in the process of agricultural intensification, attention should be paid to ecological security, especially in an oasis area with a fragile ecological environment in an arid region.