Modeling the Water and Nitrogen Management Practices in Paddy Fields with HYDRUS-1D

Abstract

Rice production involves abundant water and fertilizer inputs and is prone to nitrogen (N) loss via surface runoff and leaching, resulting in agricultural diffuse pollution. Based on a two-season paddy field experiment in Jiangsu Province, China, field water and N dynamics and their balances were determined with the well-calibrated HYDRUS-1D model. Then, scenarios of different controlled drainage and N fertilizer applications were simulated using the HYDRUS-1D model to analyze the features and factors of N loss from paddy fields. Evapotranspiration and deep percolation were the two dominant losses of total water input over the two seasons, with an average loss of 50.9% and 38.8%, respectively. Additionally, gaseous loss of N from the whole soil column accounted for more than half of total N input on average, i.e., ammonia volatilization (17.5% on average for two seasons) and denitrification (39.7%), while the N uptake by rice accounted for 37.1% on average. The ratio of N loss via surface runoff to total N input exceeded 20% when the N fertilizer rate reached 300 kg ha⁻¹. More and longer rainwater storage in rice fields under controlled drainage reduced surface runoff losses but increased the risk of groundwater contamination by N leaching. Therefore, compared with raising the maximum ponding rainwater depth for controlled drainage, optimizing N fertilizer inputs may be more beneficial for controlling agricultural diffuse pollution by reducing N loss via surface runoff and leaching. The HYDRUS-1D model provides an approach for the quantitative decision-making process of sustainable agricultural water and N management.

1. Introduction

Appropriate water management practices are indispensable in agricultural production systems to create suitable soil moisture and ensure the supply of nutrients such as nitrogen (N). Different from the water scarcity in arid areas, humid areas face the great challenge of agricultural diffuse pollution, especially in rice paddies [1]. Rice feeds more than 50% of the world’s population [2], and its high and stable yield is inseparable from the high application rate of N fertilizer.

Excessive N application is easily lost to the surrounding water environment as a result of surface runoff, which has become the major inducement of lake eutrophication and groundwater contamination [3,4,5]. Due to differences in rainfall, soil, and field management strategies in different regions, N runoff loss accounts for 5.6% to 68.2% of the total N input [6,7]. The low recovery rate of fertilizer for crops leads to significant leaching of residual N into groundwater due to irrigation or rain [8]. In particular, nitrate contributes more to groundwater pollution than any other form of N due to its fluidity [9]. In China, it was reported that more than 50% of paddies required reduced N application to meet the environmental optimum of N-fertilizer, i.e., 169–199 kg N ha⁻¹ [10]. Many studies have previously revealed the impact of reduced N application rates on crop growth or yield [11,12,13]. However, there is still a lack of effective quantitative research to assess the impact of optimizing N application on reducing agricultural diffuse pollution in paddy fields, which is unlikely to convince farmers and policymakers to pay attention to agricultural water management and optimize N fertilizer investments.

In recent years, extreme hydrological events have occurred more frequently due to climate change, especially threatening the security of agriculture production [14,15]. For rice, the introduction of alternate wetting and drying (AWD) regime, one of the widely promoted water-saving irrigation methods, regulates the field surface’s wetting or drying state and improves water productivity in rice production [16]. In humid areas, rainfall is abundant, and most crops depend on rain-fed irrigation. The uncertainty of time, frequency, and volume of rainfall can affect the movement of soil water [17] and change the transport and transformation process of N correspondingly [7]. Moreover, the paddy field experiences extreme soil moisture fluctuations under AWD conditions, which significantly affects water storage capacity and runoff output [18].

Farmers often perform artificial drainage to avoid continuous flooding of paddy fields, resulting in more than 23.9% loss of N runoff [18,19], which has been recognized as a critical source of diffuse agricultural pollution [20]. However, despite the high frequency of surface runoff in paddy fields, the characteristics of nutrient loss under controlled drainage conditions have attracted less attention. Within the rapid degradation stages (less than 10 days) of N fertilizer such as urea, it is reported that the N concentration in paddy water decreases exponentially [21,22], and significant variations of N content were observed in intra- and inter- runoff events [17]. Consequently, it is important to investigate the response relation between controlled drainage and N losses in paddy fields after rainfall in humid areas.

On-site monitoring over years is primary and helpful, but time-consuming and laborious to accurately quantify N loss loads. In particular, when analyzing, evaluating, and predicting the combination of various water and N management strategies, model simulation methods cannot be neglected and are of great significance to agro-hydrology research [23]. Several models, such as CERES-Rice [24], APSIM [25], and ORYZA [26], have been applied in rice production systems. These models can simulate the crop growth process very well, but the simulation of N transport and transformation is limited. In this regard, HYDRUS-1D has been applied extensively to simulate water, heat, and solute transport in various environments [27,28]. The sequential first-order decay chain is used to describe N transformation among urea, ammonia, and nitrate. It has the flexibility to consider the absorption of nutrients by roots and simulates the dynamics of soil water and solutes under different boundary conditions [29].

In southern China with a humid climate, frequent heavy rainfall events significantly increase surface runoff and leaching loss of solutes. Surface runoff happens particularly when rainfall exceeds paddy storage capacity, and it is determined by the height difference between the drainage outlet height (h_max, i.e., the maximum ponding rainwater depth) and the ponding water depth (h_p) [30], which are relied on manual field measurements. Agro-hydrology models such as HYDRUS-1D may be helpful in quantifying field runoff process and water and solute balances under the background of increasing frequency, intensity, and duration of extreme weather events. However, the relevant evidence is relatively scarce in this regard.

In the present study, the environmental benefits of controlled drainage and N fertilizer reduction in paddy fields were investigated in a humid climate, based on a two-year experimental observation and a well-calibrated HYDRUS model simulation. We designed different combinations of rainwater storage and N application rates and evaluated the risk of N loss via leaching and surface runoff. Under different scenarios, modeling techniques were utilized to accurately determine the elements of water and N balance in paddy fields. Research results are expected to guide water and N fertilizer management in paddy fields for sustainable agriculture.

2. Materials and Methods

2.1. Study Area and Experiment Design

The experiment was conducted at the Jiangning Water-saving Park of Hohai University, Jiangsu province, China (31°55′ N, 118°46′ E). This region has a subtropical monsoon climate. The annual average temperature, sunshine hours, and evaporation are 15.7 °C, 2017.2 h, and 900 mm, respectively. The annual average rainfall is about 1072.9 mm with an average annual 120 rainfall days, of which the period between May and September accounts for 60% of the annual rainfall. With the temporarily uneven distribution of rainfall, surface runoff from paddy fields may occur from time to time. There are 32 non-weighing drainage lysimeters in the Water-saving Park, 2.5 m × 2.0 m in length and width and 2.0 m in depth, which are separated by concrete and then sealed with waterproof coatings [31]. The plots were cultivated with rice–wheat rotation for many years and the basic properties of the test soils are listed in

Table 1. Basic soil properties of different layers.

Soil Layer (cm)	Soil Particle Fraction (%)			Bulk Density (g cm−3)	Organic Matter (%)	pH (H2O)
	Sand	Silt	Clay
0–20	40.21 (9.06)	38.22 (6.43)	21.57 (3.26)	1.36 (0.23)	1.24 (0.06)	6.94 (0.06)
20–40	39.12 (6.31)	39.16 (4.71)	21.72 (2.63)	1.40 (0.19)	1.35 (0.06)	6.97 (0.05)
40–60	38.87 (5.46)	39.86 (4.06)	21.27 (2.83)	1.43 (0.20)	1.20 (0.08)	6.85 (0.07)
60–160	40.25 (5.02)	38.12 (3.72)	21.63 (2.41)	1.48 (0.21)	-	6.80 (0.04)

Note: The standard deviations calculated from six soil samples are given in parentheses.

.

The experiments with the rice cultivar Nanjing 9108 were carried out in 2020 and 2021. Rice seedlings at 35 days old were transplanted into the lysimeters on 17 June 2020 and 4 July 2021 at 20 cm × 15 cm spacing. The N application rate was 225 kg N ha⁻¹, and urea was applied as the base fertilizer (17 June 2020 and 5 July 2021), tiller initiation fertilizer (30 June 2020 and 7 July 2021), and spikelet-developing fertilizer (2 August 2020 and 14 August 2021) at a ratio of 4:3:3. All plots were managed with practices commonly adopted by local farmers, including P and K fertilization, pesticide application, and water regimes. Specifically, water management in paddy fields followed an improved AWD technology adapted to southern China [32], and its irrigation thresholds are shown in

Table 2. Control thresholds of field water levels during the rice-growing season.

Irrigation Criteria (cm)	Returning Green	Tillering	Jointing-Booting	Heading-Flowering	Milky Ripening	Yellow Ripening
Lower threshold	1.0	−20.0	−30.0	−20.0	−30.0	Naturally drying
Upper threshold	3.0	3.0	3.0	3.0	3.0

. When the field water level dropped to the lower threshold of irrigation, paddy was irrigated with an electromagnetic flowmeter until the upper threshold was reached. On the other hand, surface drainage occurred when rainfall caused h_p to exceed drainage outlet height at 10 cm.

2.2. Measurements and Analysis

Daily meteorological data, such as air temperature, wind speed, radiation, and humidity, were measured by an automatic weather station (31°55′ N, 118°46′ E, altitude of 15.0 m) in the Water-saving Park. According to the irrigation instructions (Table 2), the field water level dynamic was measured at 9 a.m. each day by a perforated PVC pipe [33]. Irrigation of each lysimeter was controlled by a solenoid valve connected to a computer system, and the irrigation water volume was recorded by the flowmeter at the water inlet. Similarly, the amount of surface runoff was measured by a flowmeter at the outlet. In general, half an hour after the irrigation or drainage event, the field water level was monitored to ensure that it reached the lower or upper threshold.

Within two weeks after each N fertilizer application, water samples were collected at 9:00 a.m. daily for the first five days, and then every three days or after irrigation/rain events. Using a small electric vacuum pump, water samples were taken by the siphon method in the ponded layer of each lysimeter, and three different positions were taken for mixing each time. Before rice transplanting, several bottom porous PVC tubes were implanted at 20, 40, and 60 cm soil depths in the lysimeters to collect leachate [34]. Electrical water pumps were also used to extract leachate from the PVC tubes. If there was no sampling on the adjacent days, the residual leachate in the PVC tubes was drained the night before. The water samples were stored quickly at −20 °C after being collected.

The N content in water samples was determined by a spectrophotometer (P4UV, MAPADA, Shanghai, China). In specific, ammonia nitrogen (NH₄⁺-N), and nitrate-nitrogen (NO₃⁻-N) were analyzed by Naismith’s reagent UV spectrophotometry, and the hydrochloric acid acidification-UV spectrophotometry, respectively. The cumulative N loss was the product of water flux of surface runoff or percolation rate and the corresponding NH₄⁺-N and NO₃⁻-N concentrations. Moreover, leaf area index (LAI) was measured every ten days with an AccuPAR LP-80 analyzer (Decagon, Pullman, WA, USA) during rice growth. Rice samples of aboveground organs were first dried at 70 °C for 48 h and then weighed and ground. Total N contents were subjected to the Kjeldahl method and digested with H₂SO₄-H₂O₂ at 260 °C [35]. Cumulative total N uptake was dry matter weight multiplied by the measured total N content.

Soil samples were collected ten days before rice transplanting. After removal of surface debris, soil samples (in triplicate) were collected from 2 quadrats along the diagonal of each plot by soil drill. In addition, undisturbed soil cores were collected using a 5 cm-diameter ring knife to determine soil bulk density, which was based on the inner diameter of the ring knife, sampling depth, and soil drying weight. The air-dried soil samples were sieved through a 2 mm sieve and then pretreated with 6% H₂O₂ after which 10% HCl was added to remove oxides and carbonates. A laser particle size analyzer (s Malvern Panalytical Ltd., Malvern, UK) was applied to determine particle-size distribution after chemical treatment and mechanical dispersion.

2.3. Model Description

Since lateral seepage was minimized in paddy plots, the HYDRUS-1D model was applied to simulate and predict one-dimensional soil water movement, mass balance, and N dynamics in rice paddies.

2.3.1. Soil Water Flow

The Richard’s equation was used to calculate one-dimensional soil water flow [29]:

$$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}\left[K\left(h\right)\left(\frac{\partial h}{\partial z}+1\right)\right]-S_r$$

(1)

where z and t denote the soil profile coordinate (cm) and time (day) respectively; θ is the volumetric water content (cm³ cm⁻³); S_r represents the sink term of root (cm³ cm⁻³ day⁻¹); K(h) is the unsaturated hydraulic conductivity (cm day⁻¹) associated with water pressure head (cm). In addition, the mathematical relationships between θ(h) and K(h) are described as follows [36]:

$$\theta \left(h\right)=\left\{\begin{array}{cc}{\theta }_r+\frac{{\theta }_s-{\theta }_r}{{\left(1+{\left|\alpha h\right|}^n\right)}^m}\hfill & (h<0)\hfill \\ {\theta }_s\hfill & (h\ge 0)\hfill \end{array}\right.$$

(2)

$$K\left(h\right)=K_s{S}_e^{\lambda }{\left[1-{\left(1-S_e^{\frac{1}{m}}\right)}^m\right]}^2$$

(3)

where θ_s and θ_r denote saturated water content and residual water content (cm³ cm⁻³), respectively; K_s means saturated hydraulic conductivity (cm day⁻¹); and l, m, and α are empirical shape factors (-); and the S_e effective saturation is calculated by:

$$S_e=\frac{\theta \left(h\right)-{\theta }_r}{{\theta }_s-{\theta }_r}$$

(4)

2.3.2. Evapotranspiration and Root Water Uptake

Reference evapotranspiration (ET₀) rates are calculated by the Penman–Monteith equation [37] with daily weather data, and the potential evapotranspiration rates (ET_c) are the product of ET₀ and the crop coefficient (K_c). Then, LAI is used to partition evaporation and transpiration [38]:

$$E_p=E{T}_C\times e^{-Kgr\times LAI}$$

(5)

where K_gr is an extension coefficient for global solar radiation, which was set to 0.3 for rice [39]. Thus, the potential transpiration is the difference between ET_c and potential evaporation. In addition, a piecewise linear function was parameterized to describe the water stress response [40], including four pressure heads associated with root water uptake. The meanings are as follows: if h > h₁ and h < h₄, then the root is assumed to stop absorbing water; and water uptake maintain at the highest state between h₂ and h₃, while it linearly decreases (or increases) between h₃ and h₄ (or h₁ and h₂).

2.3.3. N Transport and Transformation

Three N forms in rice paddies including urea, NH₄⁺-N, and NO₃⁻-N are described as first-order chain reactions, calculated by following the convection–diffusion equation [29]:

$$\frac{\partial \theta C_1}{\partial t}=\frac{\partial }{\partial z}\left(\theta D_1^w\frac{\partial C_1}{\partial z}\right)-\frac{\partial q{C}_1}{\partial z}-{\mu }_{w,1}^{\prime }\theta C_1$$

(6)

$$\begin{gathered} \frac{\partial \theta C_2}{\partial t}+\frac{\partial \rho S_2}{\partial t}+\frac{\partial a_v{G}_2}{\partial t}=\frac{\partial }{\partial z}\left(\theta D_2^w\frac{\partial C_2}{\partial z}\right)+\frac{\partial }{\partial z}\left(a_v{D}_2^g\frac{\partial G_2}{\partial z}\right)-\frac{\partial q{C}_2}{\partial z} \\ -{\mu }_{w,2}^{\prime }\theta C_2+{\mu }_{w,1}^{\prime }\theta C_1+{\gamma }_{s,2}\rho -r_{a,2} \end{gathered}$$

(7)

$$\frac{\partial \theta C_3}{\partial t}=\frac{\partial }{\partial z}\left(\theta D_3^w\frac{\partial C_3}{\partial z}\right)-\frac{\partial q{C}_3}{\partial z}-{\mu }_{w,3}^{\prime }\theta C_3+{\mu }_{w,2}^{\prime }\theta C_2-r_{a,3}$$

(8)

where the subscripts 1, 2, and 3 represent urea, NH₄⁺-N, and NO₃⁻-N, respectively; C, S, and G denote the N concentration in liquid (mg L⁻¹), solid (mg g⁻¹), and gas (mg g⁻¹) phases, respectively; ρ is soil dry bulk density (g cm⁻³); D^w and D^g represent N dispersion coefficients for liquid and gas phase, respectively (cm² day⁻¹); q is the volumetric water flux (cm day⁻¹); µ_w is the first-order rate constant for N in the liquid phase (day⁻¹), µ′_w is the first-order rate constant for chain reaction between urea, NH₄⁺-N and NO₃⁻-N (day⁻¹); γ_s is a zero-order rate constant in the solid phase (day⁻¹); r_a is a passive root N uptake.

2.4. Initial and Time-Variable Boundary Conditions

The HYDRUS-1D model was simulated on a daily time scale starting from transplanting. As illustrated in Figure 1, the paddy soil domain with a depth of 160 cm represented the lysimeter soil column, which was divided into four layers. Since the paddy field allowed water to pond on the surface, the upper boundary condition was set as the atmospheric boundary condition with the surface layer, of which the maximum thickness (i.e., h_max) was specified as 10 cm before surface runoff was initiated. The bottom boundary condition was defined as the deep drainage boundary based on percolation flux observed at the bottom of the lysimeters. The initial condition for water flow was assigned as the shallow ponded water depth during transplanting for each season.

Initial ammonia and nitrate N concentrations in the soil profile were obtained by linear interpolation of observed data in prescriptive depths, while the initial urea concentration was the N application rate of base fertilizer divided by h_p at transplanting. Similarly, the concentration of urea taken as the upper boundary condition was calculated by the top-dressed fertilizer rate and incoming water flux at the tillering and heading stages, while the Dirichlet boundary was set as the bottom boundary condition. Since irrigation and drainage were no longer managed in the paddy field after entering the yellow-ripening stage, the simulation was started the day after transplanting until the end of the milky-ripening stage. The discretization scheme uses 111 nodes distributed in the soil profile, i.e., a spatial step of 1 cm for the upper three soil layers (0–20 cm, 20–40 cm, and 40–60 cm), and a spatial step of 2 cm for the bottom 60–160 cm soil layer.

2.5. Model Input Parameters

The van Genuchten soil hydraulic parameters were estimated by the Rosetta-Lite module and then calibrated according to the observed ponding water depth (

Table 3. The calibrated soil hydraulic parameters for different soil layers.

Depth (cm)	θr (cm3 cm−3)	θs (cm3 cm−3)	α (cm−1)	n	Ks (cm day−1)
0–20	0.066	0.418	0.0095	1.526	15.19
20–40	0.065	0.408	0.0095	1.521	7.48
40–60	0.063	0.399	0.0096	1.512	11.27
60–160	0.062	0.391	0.0108	1.484	9.08

Note: θ_r is the residual water content; θ_s is the saturated water content; α and n are the van Genuchten shape parameters; and K_s is the saturated hydraulic conductivity.

). Rosetta considered the measured particle size distribution and bulk density dataset as shown in Table 1. The Feddes model parameters, including h₁, h₂, h_3high, h_3low, and h₄, have been optimized for rice at 100 cm, 55 cm, −160 cm, −250 cm, and −15,000 cm [41], respectively. For the analysis of solute transport, the longitudinal dispersivity was initially considered equal to one-tenth of the profile depth and then calibrated [42], the molecular diffusion coefficients in free water (D_w) for NH₄⁺-N and NO₃⁻-N were 1.52 and 1.64 cm² day⁻¹, and D_g in soil air for NH₃ was 18,057.6 cm² day⁻¹. N uptake rates were calculated using linear interpolation for the observed total N content of rice between two sample time intervals. It was assumed that mineral N (NH₄⁺-N and NO₃⁻-N) could be unlimited passively absorbed by rice roots (c_max = 1000 mg L⁻¹) while urea could not (c_max = 0 mg L⁻¹). Besides, the active uptake was considered only for NH₄⁺-N and its Michaelis-Menten constant was set as 0.31 mg L⁻¹ [43].

The N first-order chain reaction processes were temperature and water content dependent, but it is a common practice to neglect these dependences [44]. The parameters of the N chain reaction were hard to measure practically, thus the parameters were initially taken from published literature [45,46,47] and then calibrated and validated (

Table 4. The calibrated parameters of N first-order chain reaction for different soil layers.

Depth (cm)	DL (cm)	Kd (cm3 g−1)	Kv (day−1)	Kn (day−1)	Kdn (day−1)	Km (day−1)	Kh (day−1)
0–20	7.5	3.5	0.03	0.08	0.25	0.0045	0.74
20–40	5.5	3.5	0.03	0.02	0.05	0.0045	0.74
40–60	6.2	3.5	0.01	0.03	0.01	0	0
60–160	16	3.5	0.01	0.02	0.01	0	0

Note: D_L is the longitudinal dispersivity; K_d is the distribution coefficient for ammonia N; K_v is the ammonia volatilization coefficient; K_n is the nitrification rate; K_dn is the denitrification rate; K_m is the comprehensive production rate of NH₄⁺-N representing mineralization/immobilization; and K_h is the hydrolysis rate.

). It was reported that mineralization/immobilization was mostly active in the root zone, while ammonia volatilization only happened on the field surface [47]. Besides, the parameters of nitrification and denitrification were initially defined as the same in the soil column and then calibrated for each layer. The NH₄⁺-N distribution coefficient (K_d) between its liquid and solid phase was set as 3.5 cm³ g⁻¹, which was consistent with previous literature by Tan et al. [34], Hanson et al. [44], and Dash et al. [48], while urea and NO₃⁻-N only existed in the liquid phase (i.e., K_d = 0 cm³ g⁻¹).

The well-fit parameters were obtained through trial and error until the simulated model results closely matched the observed ponding water depth and N concentrations. Three statistical indicators including determination coefficient (R²), Nash–Sutcliffe modeling efficiency (NSE), and root mean square error (RMSE) were calculated and served as a measure of the agreement between observed and simulated data (

Table 5. Model calibration and validation results of h_p and NH₄⁺-N and NO₃⁻-N concentrations in different soil depths. h_p denotes ponding water depth of paddy field.

Period	Type (Units for RMSE)	Depth (cm)	RMSE	NSE	R2
Calibration 2020	hp (cm)	0	0.382	0.980	0.981
	NH4+-N concentration (mg L−1)	0	4.608	0.451	0.669
		20	0.729	0.609	0.712
		40	0.399	0.729	0.845
		60	0.002	0.605	0.608
	NO3−-N concentration (mg L−1)	0	2.310	0.756	0.895
		20	1.139	0.911	0.921
		40	0.461	0.880	0.881
		60	0.230	0.947	0.947
Validation 2021	hp (cm)	0	0.373	0.981	0.983
	NH4+-N concentration (mg L−1)	0	5.570	0.393	0.537
		20	0.453	0.576	0.704
		40	0.436	0.773	0.777
		60	0.011	0.440	0.517
	NO3−-N concentration (mg L−1)	0	2.568	0.767	0.881
		20	1.407	0.892	0.893
		40	0.190	0.960	0.962
		60	0.116	0.952	0.957

).

2.6. Model Simulation and Evaluation

• Scenarios Design and Model Runs

According to the existing drainage outlet depth (10 cm), the maximum ponding rainwater depth (h_max) was changed to 5–15 cm, in steps of 1 cm. N fertilizer was applied at eleven levels in the range of 100 to 350 kg ha⁻¹ increased by 25 kg ha⁻¹. In particular, N fertilizer was applied at 40%, 30%, and 30% at transplanting, tillering, and panicle initiation stages, respectively. For each scenario of combined N application and controlled drainage, the simulation was for two typical years (2020 and 2021 rice-growing seasons). Therefore, there were 11 × 11 = 121 treatments and 121 × 2 = 242 model runs in total. The coupling environmental effects between different N fertilizer levels and controlled drainage thresholds were analyzed according to the N loss load in the percolation from the root zone and surface runoff under different scenarios calculated by the HYDRUS-1D model.

• Model Performance Criteria

$$RMSE=\sqrt{\frac{1}{n}\sum _1^n{\left(P_i-O_i\right)}^2}$$

(9)

$$NSE=1-\frac{\sum {\left(P_i-O_i\right)}^2}{\sum {\left(O_i-\overline{O}\right)}^2}$$

(10)

$$R^2={\left\{\frac{{\displaystyle \sum _{i=1}^n}\left(O_i-\overline{O}\right)\left(P_i-\overline{P}\right)}{\sqrt{{\displaystyle \sum _{i=1}^n}{\left(O_i-\overline{O}\right)}^2}\sqrt{{\displaystyle \sum _{i=1}^n}{\left(P_i-\overline{P}\right)}^2}}\right\}}^2$$

(11)

where n denotes the number of results and subscript number i denotes its order; P_i and O_i denote the simulated and observed values, respectively, with the bar above letters indicating their means.

3. Results and Discussion

3.1. Model Calibration and Validation

The observed data during the 0–90 DAT in the 2020 season were used for model calibration, and then the HYDRUS-1D model with calibrated parameters was validated during the 0–80 DAT in the 2021 season. The soil hydraulic parameters were estimated by the Rosetta approach and then calibrated and validated according to the observed values of h_p, and then the calibrated water model was applied to determine the nitrification and denitrification rates at 20, 40, and 60 cm depths. The agreement between observed and simulated h_p during calibration was good (Table 5), with R² = 0.981, NSE = 0.980, and RMSE = 0.382 cm. And the validation performance of calibrated water model was similarly good, with R² = 0.983, NSE = 0.981, and RMSE = 0.373 cm. The HYDRUS-1D model does not consider the preferential flow due to macropores and cracks under AWD irrigation, which would cause an increase in percolation and a decrease in surface drainage to a certain extent, especially during the periods of dry and wet alternation of paddy fields [34], which may reduce the accuracy of simulating soil water flow movement.

For NH₄⁺-N, the model performed better for subsurface concentrations than for the surface (Table 5), with R² = 0.537–0.669 (NSE = 0.393–0.451 and RMSE = 4.608–5.570 mg L⁻¹) for surface concentration, and R² = 0.517–0.845 (NSE = 0.440–0.773 and RMSE = 0.002–0.729 mg L⁻¹) for subsurface concentration. The observed NH₄⁺-N concentrations at 0 cm and 60 cm had relatively low R² with respect to the simulation, but the simulated trend in the soil column was consistent with observation. In contrast, the model performance on NO₃⁻-N was generally better; the simulation accuracy of NO₃⁻-N was higher than that of NH₄⁺-N, and NO₃⁻-N concentration values decreased from the surface along the soil profile (Figure 2). Moreover, the simulation accuracy was better with the increased soil depth, with R² = 0.881–0.895, NSE = 0.756–0.767, and RMSE = 2.310–2.568 mg L⁻¹ for surface concentration, and R² = 0.881–0.962, NSE = 0.880–0.960, and RMSE = 0.116–1.407 mg L⁻¹ for subsurface concentration. Overall, the present study showed that the observed and simulated values of NH₄⁺-N or NO₃⁻-N concentrations were in good agreement compared with similar studies in rice paddies [46,48].

As illustrated in Figure 2a, the values of simulated NH₄⁺-N concentration were slightly underestimated with a slope of 0.568 (R² = 0.783). It is often difficult to improve the accuracy of simulated N concentrations because paddy field conditions such as pH, soil aeration, temperature, and C/N ratio may have individual or combined effects on N transformation [49]. When urea was hydrolyzed to NH₄⁺-N, the accompanying volatilization of NH₃ that was easily soluble in water was also determined by Naismith’s reagent UV spectrophotometry, which may make the observed surface concentration of NH₄⁺-N greater than the simulated value. Moreover, at 60 cm depth, it may be that the sampling pipe has been punctured through the plow pan so that the observed NH₄⁺-N concentration was larger and more fluctuating than the simulated value. With a slope of 1.078 (R² = 0.914), simulated NO₃⁻-N concentration matched measured values well (Figure 2b). Although previous studies have confirmed that it was reliable for HYDRUS-1D to be applied in lowland rice fields, some problems still existed including a sampler, macropores, earthworm burrows, roots, and soil cracks [50,51], which would affect the accuracy of water and N transport calculations.

3.2. Analysis of Field Water Dynamics and Balance

Figure 3a,b shows the dynamics of precipitation, irrigation, surface runoff, and ponded depth during the 2020 and 2021 rice-growing seasons, respectively. Apparently, the overall trend of observed h_p was in good agreement with simulated values, but differences between simulated and observed values during 2020 and 2021 mainly occurred during rainstorm events and alternate drying and wetting periods. During the 0–35 DAT in 2020, frequent heavy rainfall lasting for about one month kept h_p at more than 4 cm. In addition, three drainage events exceeding the upper limit of paddy rainwater storage resulted in a total surface runoff of 6.5 cm. Due to the large evapotranspiration of paddy fields on continuous sunny days, six irrigation events in the periodic dry/wet cycles occurred in the latter half of the season. In contrast, in the 2021 season, a rainstorm of 14.19 cm per day on the 24 DAT raised h_p to a maximum depth of 10 cm, causing surface runoff which reached 7.9 cm within one day. Under the regime of controlled drainage, the paddy field experienced long-term flooding and anaerobic conditions (Figure 3). The alternating occurrence of drought and flood may have an adverse impact on the rice growth process [52]. How to optimally control the rainwater storage amount and time to reduce the impact on rice yield under water-saving irrigation paddies merits further research.

As shown in Figure 4, the cumulative precipitation in the 2021 season was 59.9 cm, which was 5.1% lower than that in the 2020 season, but the rainwater use efficiency (i.e., surface runoff/rainfall) in the 2021 season was 15.4% higher than that in the 2020 season with irrigation reduced by 35.4%. According to simulated values, h_p above 3 cm during the 2020 and 2021 seasons were 47 and 43 days, respectively, and the corresponding deep percolation was 36.2 and 25.8 cm, respectively. Moreover, simulated ET (evaporation + transpiration) in 2020 and 2021 accounted for 49.3% and 52.6% of total water input (TWI) respectively, and deep percolation accounted for about 42.1% and 35.5% of TWI, respectively. Frequent heavy rain and more irrigation events during 2020 resulted in much more percolation water. Actually, the paddy plow pan with relatively low buffers the vertical soil water movement during the dry/wet alternating stage [53], making it possible to couple water-saving irrigation and controlled drainage in paddy fields. The water balance errors were 0.8% and 1.7% of TWI, indicating that it was reliable in applying HYDRUS-1D to analyze the water balance in rice paddies. With less rainwater loss via surface runoff, controlled drainage applied in paddy fields could utilize freshwater resources sufficiently for irrigation in the humid climate.

3.3. Analysis of N Concentration Dynamics and N Balance

Figure 5 shows the dynamics of surface NH₄⁺-N and NO₃⁻-N concentrations for two seasons, and the overall trend of the simulation was in good agreement with the measured value. In particular, the N concentration increased abruptly after the soil surface dried up and encountered irrigation or rainfall, whose peak values were difficult to capture by field sampling and monitoring. Frequent intense rainfall also reduced the accuracy of the model simulation results matching the observed N concentrations. The distribution of N leaching at the depth of 0–100 cm soil profile is shown in Figure 6. N leaching occurred mainly within two weeks after fertilizer application, especially under conditions of large rainwater storage at the paddy field surface. In both seasons, two peak values of N concentrations were observed after urea fertilization, tillering and spikelet-developing fertilizer, respectively. The concentrations of NH₄⁺-N and NO₃⁻-N at the upper 0–10 cm soil were higher than 3.0 and 9.0 mg L⁻¹ within 1–2 weeks after fertilization, respectively. The NO₃⁻-N values throughout the soil column were higher than that of NH₄⁺-N, and the NO₃⁻-N leaching was more frequent than that of NH₄⁺-N. In addition, the plow pan with low permeability at depths of 20–40 cm prevented nutrients from leaching deeper soil layers (Figure 6), which was fundamental to rice water-saving irrigation.

N dynamics in paddy fields are uncertain and time-varying affected by the wetting and drying states of soil [17]. Urea was hydrolyzed to NH₄⁺-N first, then diffused mainly to the upper 0–20 cm soil and continuously provided an N source for nitrification. In the latter half of each season, a concentration decline of NH₄⁺-N was observed under aerobic conditions of alternating dry and wet in paddy fields. It should be noted that compared with NO₃⁻-N, NH₄⁺-N was prone to adsorbing to negatively charged soil, thus significantly retarding its leaching [54]. Long-term rainwater storage may reduce oxygen in the rhizosphere, which prevents nitrification of NH₄⁺-N, as well as its leaching [55]. Moreover, relatively high NO₃⁻-N leaching continued before the soil wetting/drying stage. The concentration of NO₃⁻-N above the plow pan remained above 3 mg L⁻¹ for two seasons (Figure 6), and the difference in surface flooding duration affected the soil depth of N loads leaching. Considering that rice roots were predominantly distributed in the upper 60 cm depth of the soil column, the transport and transformation process of N in the upper soil layer affected the efficiency of nutrient uptake by rice.

Quantitative and precise analysis of N balance is beneficial to understanding N fertilizer use efficiency and loss pathways. As illustrated in Figure 4, during the 2020 and 2021 seasons, N gaseous loss from the simulated 160 cm soil column accounted for more than half of total N input (TNI) on average, i.e., ammonia volatilization (17.5%) and denitrification (39.7%), while the N absorbed and utilized by crops accounted for 37.1% of TNI on average. Except that N runoff loss was 5.8 and 16.0 kg ha⁻¹ for the 2020 and 2021 seasons respectively, other N output components were close. Both N leaching and soil N depletion accounted for less than 5% of TNI on average. The total N balance errors of simulations were less than 3.8% of TNI for the 2020 and 2021 seasons. Typically, nitrification and denitrification are two important processes in paddy fields, and ammonia volatilization and N₂O emission are important ways of N loss in the paddy field [56,57,58]. The N balance results were similar to previous studies [47,49], which were affected by the degree of soil saturation, oxygen, and temperature [59]. Additionally, rice roots and microorganisms, variably distributed in the soil, significantly affect the N transformations, which makes the N distribution in the soil profile more complex. In terms of the N balance system of the entire rice field, TNI component such as wet deposition was neglected in the calculation.

3.4. Scenarios of Paddy Controlled Drainage and N Fertilizer Application

Generally, paddy fields have the characteristics of bunds and plow pans. The former is established to allow surface water ponded to maintain the humid environment for lowland paddy, and the latter can retard the water and nutrients leaching from the root zone due to its low hydraulic conductivity. Figure 7 shows surface runoff and percolation from the root zone at a depth of 60 cm under different thresholds of h_max for controlled drainage. The percolation of the root zone greatly increased when h_max ranged from 5 to 10 cm, and then the increment of percolation slowed down when h_max continued to increase. In the seasons of 2020 and 2021, when h_max was increased from 5 cm to 10 cm, surface runoff loss of paddy field was reduced by 45.3% and 56.7%, respectively, and the corresponding leakage was increased by 17.0% and 37.8%, respectively; when h_max was increased from 10 cm to 15 cm, surface runoff loss of paddy field was reduced by 69.0% and 63.2%, respectively, and the corresponding leakage was increased by 6.5% and 5.0%, respectively. Over the two typical years, the maximum surface runoff and percolation reached 17.8 cm and 41.7 cm, respectively. Downward leaching from the root zone represented the largest portion of the available water for rice. Based on mini-tensiometers and double-ring infiltrometer measurements [60], it was reported that the infiltration rate of the paddy field increased by 1.5 times when h_p increased from 6 cm to 16 cm. Leaching predominantly depended on the soil hydraulic conductivity of each soil layer and the overall pressure head gradient that changed with surface ponding depth [53].

Runoff and leaching are the direct pathways of N load loss from farmland to the adjacent water environment [61]. Figure 8 shows that the changes in cumulative N leaching loss and surface runoff loss are similar under all scenarios of controlled drainage and N fertilizer application. N loss via surface runoff was very small and showed little variation when h_max was greater than 12 cm in the 2020 season and 8 cm in 2021. Under the current N application rate (225 kg N ha⁻¹), when h_max was set to 15 cm, the average N leaching loss was 11.10 kg ha⁻¹ and the average N runoff loss was 0.44 kg ha⁻¹ for two typical years. When h_max was lower than 8 cm, surface runoff was more sensitive to the threshold of controlling drainage. The proportion of N runoff loss to TNI exceeded 20% when the N application rate reached 300 kg ha⁻¹. Meanwhile, the risk of N leaching loss was relatively low and stable when the strategy of low fertilizer rate and low rainwater storage was adopted. In both the 2020 and 2021 seasons, N leaching loss reached above 12 kg ha⁻¹ when the N application rate was between 250 and 350 kg ha⁻¹ and h_max was higher than 10 cm. N loss via leaching and surface runoff showed a positive trend with increasing N application levels. Thus, our study suggested that it was difficult for controlled drainage to reduce N loss via surface runoff and leaching at the same time while controlling the N application rate was more sensitive to reducing the non-point pollutants of N.

Unreasonable agricultural water management practices, such as uncontrolled over-irrigation and leakage loss from irrigation channels, may make the lowland paddies prone to flooding and waterlogging [62]. Although a large threshold of h_max can be set in the rice water-saving irrigation regime to reduce surface runoff, more and longer rainwater storage in paddy fields increased the loss of deep percolation, which may increase the risk of N leaching. The higher soil surface water pressure was accompanied by more vertical water flow, as well as N leaching. Moreover, long flooding and anaerobic conditions may have a negative impact on rice growth and yield [52]. However, HYDRUS-1D faced difficulties in assessing the impact of flooding on crop production due to its limited number and accuracy of modules. In further investigation, it is necessary to couple hydrologic and crop growth models to analyze the comprehensive processes of soil water movement, nonpoint nutrient loss, and crop growth, and develop a more robust water-saving and environmentally friendly rice production technology.

4. Conclusions

Field water dynamics and N transport and transformation of a rice paddy have been monitored and simulated with the HYDRUS-1D model. The calibrated results showed that the model has better accuracy in simulating h_p and NO₃⁻-N concentration, followed by the simulation of NH₄⁺-N. Evapotranspiration and deep percolation were two dominant losses of TWI across two seasons, while ammonia volatilization and denitrification from the entire soil column accounted for more than half of TNI. Plow pan with low permeability reduced the downward migration of water and solute, especially NH₄⁺-N. The proportion of N runoff loss to TNI exceeded 20% when the N application rate reached 300 kg ha⁻¹. Additionally, it should be noted that more and longer rainwater storage in paddy fields increased the loss of deep percolation, which increased the risk of N leaching. Therefore, it was difficult to achieve optimal control of N loss via both surface runoff and leaching at the same time under controlled drainage, and reducing the N fertilizer level was more effective than increasing the allowable rainwater storage threshold of controlled drainage. It merited further investigation on the comprehensive effects of controlled drainage and N application on rice production and diffuse agricultural pollution with the coupled hydrologic and crop growth models.