HYDRUS-2D simulations of nitrate nitrogen and potassium transport characteristics under fertilizer solution infiltration of furrow irrigation Open Access

Understanding the characteristics of soil solute transport is fundamental to the design and management of furrow irrigation systems. However, the study of this problem is complicated because of the nonlinear relationship between soil capillary force and gravity at different locations during furrow irrigation, which can be affected by soil physical properties and different irrigation technical factors, and the water potential gradient of fertilizer solution infiltration varies at different locations. Therefore, some researchers (Li et al. 2004; Ebrahimian et al. 2012b) have analyzed the effects of different factors on the transport characteristics of water content and fertilizer (such as nitrate nitrogen NO₃⁻-N and ammonium nitrogen NH4⁺-N) in soil wetting patterns based on fertilizer solution infiltration, and numerous empirical models have been established, which can be used to estimate cumulative infiltration, water content of wetting patterns, and NO₃⁻-N and NH4⁺-N distributions. Due to the influences of soil texture and initial condition, the distributions of different types of ions in the soil exhibit obvious differences. Silberbush & Barber (1983) indicated that potassium ions (K⁺) have strong adsorption characteristics in soil. Hanson et al. (2006) simulated and analyzed the transport processes of K⁺ under drip irrigation, and the results proved that K⁺ can easily accumulate on soil surfaces. Liu et al. (2017) studied results that showed that the distribution of NO₃⁻-N was similar to that of soil water content, whereas K⁺ was concentrated around the Moistube. Experimental research methods can reveal the transport characteristics of soil water and fertilizer, but experimentation is time-consuming, and most research results are limited because they are obtained under specific experimental conditions.

Researchers use various numerical simulation methods to study soil water and fertilizer transport characteristics. Zerihun et al. (2005a, 2005b) established a mathematical model of soil water flow and solute transport for border fertigation and verified that the model accurately simulated the soil water content and solute transport. Doltra & Muñoz (2010) simulated NO₃⁻-N transport processes in soil under drip fertigation, and the correlation coefficient between their simulated and measured NO₃⁻-N content was higher than 0.76. Tan et al. (2015) simulated the transport characteristics of soil water, NO₃⁻-N, and NH₄⁺-N in paddy fields, and the results showed that the simulated values were consistent with the measured values. Li et al. (2019) simulated soil water and nitrogen transport processes under drip irrigation, and the relative error between the simulated and measured soil water content was within 10%; the errors for NO₃⁻-N and NH₄⁺-N were within 20%. The aforementioned studies proved that numerical simulation can be used to study the characteristics of fertilizer solution infiltration with high accuracy. Although many studies have been done on solute transport characteristics during fertilizer solution infiltration (Deb et al. 2015; Iqbal et al. 2016; Brunetti et al. 2018; Ranjbar et al. 2019), few publications have considered the effect of the changes in fertilizer concentration on soil hydraulic parameters, and analyzed the influences of various factors (such as initial soil water content, fertilizer concentration, and water depth) on the solute transport characteristics.

Therefore, the objectives of this study were: (1) to determine the soil hydraulic and solute transport parameters by inverse solution with HYDRUS-2D from infiltration data of different potassium nitrate (KNO₃) concentrations under furrow irrigation, and then verify the reliability of the parameters; and (2) to analyze the effect of initial soil water content (θ₀), KNO₃ concentration, and water depth (h₀) on the NO₃⁻-N and K⁺ transport characteristics by HYDRUS-2D.

Soil samples were taken from Yangling District (107°55′50″–108°07′50″ E and 34°14′30″–34°19′00″ N), Shaanxi Province. The terrain of Yangling District is divided into three terraces from south to north. The altitude of the first terrace is 420–430 m, and the soil texture is sandy loam. The altitudes of the second and third terraces are 450–485 m and 515–540 m, respectively, and the soil texture is clay loam. According to the topographical features of Yangling District, soil samples were collected from typical fields cultivated all year round in the first and third terraces at a depth of 0–60 cm (the depth where the roots of crops such as wheat and corn are mainly distributed). The residual water content of air-dried soil samples was measured using a weighing method. The soil texture was determined by a Mastersizer 2000 particle size analyzer. The contents of NO₃⁻-N and K⁺ in the soil samples were determined by an ultraviolet spectrophotometer and a flame photometer, respectively. The results are listed in Table 1.

KNO₃ is a potassium–nitrogen compound fertilizer with stable physical and chemical properties. It can supplement the potassium and nitrogen required for crop growth at the same time, and it is widely used in crops and fruit trees. Therefore, KNO₃ was selected as the experimental fertilizer in this study. It was dissolved in water to carry out fertilizer solution infiltration experiments. The soil bulk density (ρ) of the first and third terraces in Yangling District varied from 1.40 to 1.50 g cm⁻³ and 1.30 to 1.40 g cm⁻³, respectively (Nie et al. 2017). Therefore, two bulk density cases were designed for sandy loam and clay loam in this study, namely 1.40 and 1.50 g cm⁻³ for sandy loam and 1.30 and 1.40 g cm⁻³ for clay loam. Based on the ranges of KNO₃ concentration set by previous studies (Zipelevish et al. 2000; Wang et al. 2010; Liu et al. 2017), the KNO₃ concentration was tested, in this study, at four gradients: 0, 250, 600, and 900 mg L⁻¹. Furthermore, θ₀ and h₀ were selected as the experimental design factors to carry out infiltration experiments of furrow irrigation (Table 2).

The air-dried soil sample was filtered through a 2 mm sieve, and a 5 cm layer of soil from the sieved material was placed in a soil box; the size of the soil box was 80 × 5 × 100 cm (length × width × height). A furrow was excavated along the side of the soil box, and that furrow had the trapezoidal section commonly used in northern China, with a maximum depth of 15 cm, bottom width of 20 cm and a 1:1 side-slope. Because the furrow was axisymmetric, only a half-section was used in the experiments (the range of ABCDO), as shown in Figure 1(a). The infiltration time and cumulative infiltration were recorded. After the experiments of infiltration, soil samples were collected from different locations (Figure 1(a)), and soil samples were divided into two parts. Some soil samples were measured for water content through a weighing method, and other soil samples were subjected to an ultraviolet spectrophotometer and flame photometer to determine their contents of NO₃⁻-N and K⁺, respectively.

Boundary conditions (Figure 1(b)): DE and EG represent the variable pressure head boundaries, and a constant water depth and solution concentration are maintained. Because of the short durations of the infiltration experiments, the scenario posits no rainfall and no plant cover on the surface, and it is considered to have no external solute supplementation or loss. Therefore, GF and FA represent the atmospheric boundaries, and evaporation can be neglected to create zero-flux boundaries. The left boundary AB is the symmetry axis, and the horizontal flux is 0. The right boundary CD is the symmetry axis of the adjacent furrow, which is regarded as a zero-flux surface. The lower boundary BC contains drainage holes. During the experiment, the wetting front does not reach the base of the soil box and does not affect the process of fertilizer solution infiltration, and the constant initial conditions are established to create a free drainage boundary.

On the basis of the initial (Equations (7) and (8)) and boundary conditions (Figure 1(b)), Equations (1) and (4) are solved by HYDRUS-2D by finite elements, which can define the computational grid automatically. The smaller finite elements are recommended in regions where rapid changes may occur of water content and pressure (such as at the infiltrating surface) and larger ones elsewhere. The resulting computational grid consisted of 2,446 nodes (not shown in Figure 1(b)).

This inverse solution procedure has been successfully applied by several researchers (Abbasi et al. 2003; Verbist et al. 2009; Ebrahimian et al. 2012a) to determine soil hydraulic and solute transport parameters. Because the θ_r has little effect on soil water movement (Wang et al. 2013), θ_r was set as the water content of air-dried soil in this study. Nitrogen in the KNO₃ solution is mainly in the form of NO₃⁻ and is negatively charged, and the adsorption of NO₃⁻ can be neglected. The soil wetting pattern has a high water-content in fertilizer solution infiltration, and the nitrification of nitrogen is suitable for aerobic environments; thus, the nitrification of NH₄⁺-N can be neglected. Furthermore, the infiltration experiments were brief (the longest time = 6 h), and consequently the mineralization and denitrification of nitrogen in the soil can be neglected in this study. The temperature during the infiltration experiments was approximately 20–25 °C. Siyal et al. (2012) reported that D_w = 0.00113 cm² min⁻¹ for NO₃⁻-N at 20–23 °C, and Ren et al. (2013) wrote that D_w = 0.00395 cm² min⁻¹ for K⁺, so the D_w of NO₃⁻-N and K⁺ can be assigned previously published values in the solute transport parameter inversion process of this study. Therefore, the inverse solution was applied to four soil hydraulic parameters, including θ_s, n, K_s, and a, and three solute transport parameters, including α_k of K⁺, D_L and D_T of NO₃⁻-N and K⁺. The inverse optimization method simultaneously uses all measured data, i.e., soil water contents and NO₃⁻-N and K⁺ concentrations.

The cumulative infiltration per unit area curve (cumulative infiltration per length divided by wetted perimeter) under different KNO₃ concentrations is shown in Figure 2. From the measured soil water contents and concentrations of NO₃⁻-N and K⁺ at the end of the infiltration experiments of Nos. T1–T8 (Table 2 and Figure 2), the soil hydraulic and solute transport parameters were determined by inverse solution with HYDRUS-2D, and the results are listed in Table 3.

The cumulative infiltration per unit area showed an increasing trend with the increase of KNO₃ concentration. For example, for the cases where KNO₃ concentrations were 250, 600, and 900 mg L⁻¹, the cumulative infiltration per unit area of clay loam and sandy loam soil textures increased by 17.05%, 28.65%, 34.71% and 12.20%, 16.19%, 20.45% compared with 0 mg L⁻¹ at the end of the infiltration experiments, respectively (Figure 2). The soil hydraulic parameter inversion results showed that the parameters of θ_s, n, and K_s increased with the increase of KNO₃ concentrations, but a showed a decreasing trend (Table 3), which is consistent with the results of Feng (2017). The reasons may be that K⁺ has strong adsorption (Silberbush & Barber 1983; Hanson et al. 2006; Ren et al. 2013), which exchanges with the ions on the surface of the soil colloid, changing the properties of the soil colloid and thus affecting the pore structure of the soil. Moreover, the KNO₃ solution has a higher solute potential than the water, and the soil–water potential difference at the interface of the wetting front is larger for the KNO₃ solution than for the water, which increases the infiltration capacity of the soil (Figure 2). This explanation would be consistent with the research conclusions of Liu et al. (2017). Meanwhile, it should be noted that the parameters of θ_s, n, and K_s were unlikely to increase continuously with the increase of KNO₃ concentrations, and a cannot decrease continuously. The reason may be that with the increase of KNO₃ concentrations, the pore structure of the soil is changed and the solute potential at the interface of the wetting front is increased, which makes the soil infiltration capacity gradually increase. However, with the continuous increase of KNO₃, the viscosity of water increases, that is, the fluidity of water becomes worse, which will reduce the soil infiltration capacity. Because only four concentrations of KNO₃ were considered in this study, and the concentrations are relatively small (the maximum concentration of KNO₃ is 900 mg L⁻¹), what leads to this process cannot be reflected upon, and further research is needed on this issue.

The solute dispersivities mainly depend on soil properties, and these parameters have scale characteristics. Beven et al. (1993) and Cote et al. (2003) wrote that D_L can be taken as 1/10 of the simulated length. The results of Šimůnek & van Genuchten (2008) and Siyal et al. (2012) proved that D_L is generally between 0.1 and 15 cm, and D_T takes 1/100–1/5 of D_L according to soil properties and can obtain a satisfactory solute transport simulation result. As isotropic soil is used in this study, and the D_L and D_T of K⁺ and NO₃⁻-N have small differences, which were determined by inverse solution with HYDRUS-2D, so the mean values of D_L and D_T were adopted under the same bulk density of each soil texture, i.e. the ranges of D_L are 4.03–4.66 cm and 1.26–1.95 cm for clay loam and sandy loam, respectively, and the value of D_T is approximately 1/20 of D_L (Table 3). These findings are basically consistent with the aforementioned results, indicating that the inversion results of D_L and D_T of NO₃⁻-N and K⁺ are reasonable in this study.

The soil adsorption process of K⁺ can be described as a linear equation (Hanson et al. 2006; Grecco et al. 2019); that is, α_k is a constant value, which is independent of the solution concentration. In addition, the adsorption of K⁺ is mainly related to soil organic matter and particle content (Long et al. 2001), and thus the α_k of K⁺ is taken as the inversion average value under the same soil texture, namely the α_k of K⁺ is 4.47 and 3.19 cm³ g⁻¹ in clay loam and sandy loam, respectively. The inversion results of α_k are close to α_k = 1.533 cm³ g⁻¹ of silt loam (Ren et al. 2013), but this result differs greatly from the 28.7 cm³ g⁻¹ (soil texture is loam) used in the research process of Hanson et al. (2006). The reasons may be that field experiments were used by Hanson et al. (2006), and numerous external factors are relevant. For example, salt may accumulate in the surface soil with water evaporation, and the solute that is absorbed by plant roots may affect the α_k of K⁺.

To verify the reliability of soil hydraulic and solute transport parameters determined by inversion solution, the parameters were input into HYDRUS-2D (Table 3) and the verification experiments of Nos. V9–V12 (Table 2) were simulated. The infiltration time used in the simulation process was consistent with the infiltration time of the corresponding experiment treatment. The simulation results and the measured values are shown in Figures 3 and 4.

The mean values of RMSE, MPBE and MBPRE of cumulative infiltration and soil water content between the measured and simulated values for the verification experiments of Nos. V9–V12 were 0.66 cm, 0.24%, 3.45%, and 0.027 cm³ cm⁻³, 1.04%, 6.04%, respectively (Table 4). The results indicated that the simulated cumulative infiltration and soil water content exhibited good consistency with measured values (Figure 3). And analysis of the causes of errors was that on the one hand that the soil hydraulic parameters were used to simulate verification experiments of Nos. V9–V12, these parameters by inverse solution from sandy loam and clay loam soil texture infiltration experiments of Nos. T1–T4, respectively (Tables 2 and 3), but some differences existed in h₀ and θ₀, resulting in errors of cumulative infiltration and soil water content between the simulated and measured experiments. On the other hand, measurement errors are inevitable during such experiments. Overall, the results indicated that the soil hydraulic parameters determined by the inversion solution are reliable and those parameters can be used to simulate the soil water movement process of furrow irrigation.

The simulated concentrations of NO₃⁻-N and K⁺ were compared with the measured values, and they agreed (Figure 4). The mean values of RMSE, MPBE and MBPRE of NO₃⁻-N and K⁺ concentrations between the measured and simulated values for the verification experiments of Nos. V9–V12 were 0.062 g cm⁻³, −0.18%, 9.34%, and 0.066 g cm⁻³, 0.05%, 10.54%, respectively (Table 4). The accuracy of the simulation results was slightly worse than the soil water movement simulation results, which was consistent with the research results of Ma et al. (2004) and Phillips (2006). The reasons may be that physicochemical reactions such as mineralization, nitrification, and denitrification of nitrogen were neglected in the process of NO₃⁻-N transport and parameter inversion in this study, and the effects of soil bulk density, water content, water depth, and fertilizer concentration on α_k of K⁺ were ignored. In addition, solute transport in soil is more complicated than water movement, and this study's description of the transport of convection dispersion equation was simplified. Considering the aforementioned reasons and the actual situation of the experiments, it can be concluded that the simulation accuracy of solute transport was reasonable; thus, the solute transport parameters determined in HYDRUS-2D are relatively reliable.

Numerical simulations were used to analyze the effect of θ₀, KNO₃ concentration, and h₀ in NO₃⁻-N and K⁺ transport characteristics by HYDRUS-2D, based on inversion solution of soil hydraulic and solute transport parameters (Table 3). Because of the similarity of the distribution characteristics of NO₃⁻-N and K⁺ for clay loam and sandy loam, the following analysis takes clay loam as an example.

The distribution of the NO₃⁻-N and K⁺ concentrations of the soil wetting pattern profile at the end of infiltration are shown in Figure 5, given the following values: KNO₃ concentration C = 250 mg L⁻¹, water depth h₀ = 10.5 cm, soil bulk density ρ = 1.30 g cm⁻³, initial water content θ₀ = 0.130, 0.156, and 0.208 cm³ cm⁻³ of clay loam.

The θ₀ influences the transport characteristics of NO₃⁻-N in the wetting pattern (Figure 5). At the end of infiltration (t = 360 min), for the cases where θ₀ was 0.156 and 0.208 cm³ cm⁻³, the vertical distribution depth of NO₃⁻-N increased by 5.28% and 13.53% compared with 0.130 cm³ cm⁻³, and the horizontal direction increased by 4.79% and 9.24%, respectively. The reasons may have been that the soil suction (mainly the matrix potential) decreased with the increase of θ₀, which may have slowed the wetting front movement. However, compared with θ_s, the water shortage (θ_s−θ₀) in the soil decreased with any increase in θ₀, which accelerated the wetting front movement. The aforementioned two reasons work together, which makes the soil wetting pattern distribution range increase with the increase of θ₀, which is consistent with the results of Lazarovitch et al. (2009) and Nie et al. (2015). Because the adsorption of NO₃⁻-N was not considered in this study, NO₃⁻-N migration was mainly affected by soil water movement; thus, the distribution range of NO₃⁻-N in soil was consistent with the distribution of soil water (wetting pattern), and the overall distribution range of NO₃⁻-N increased with the increase of θ₀. Furthermore, the results also proved that the high value region of NO₃⁻-N concentration had a downward trend with any increase of θ₀. For example, the high value regions of NO₃⁻-N concentration under θ₀ were 0.156 and 0.208 cm³ cm⁻³, moved downward by approximately 3.8 and 6.6 cm compared with 0.130 cm³ cm⁻³, respectively; this further confirmed that the characteristics of NO₃⁻-N easily move with soil water. Therefore, the results suggested that in the range of soil water content suitable for crop growth, fertigation should be avoided when the soil water content is high to minimize the risk that NO₃⁻-N might leach deep into the soil.

K⁺ mainly accumulates near the bottom of furrows (Figure 5), which indicates that θ₀ has little effect on the distribution range of K⁺. In this study, the distribution depths of K⁺ are mainly between 23 and 26 cm in the vertical direction, and within 25 and 27 cm in the horizontal direction under the conditions of θ₀ ranging from 0.130 to 0.208 cm³ cm⁻³, respectively. The reason is that the soil wetting pattern movement is accelerated with the increase of θ₀, which shortens the time of colloid adsorption of K⁺ in the soil, so the K⁺ distribution range is slightly increased.

The distributions of the NO₃⁻-N and K⁺ concentrations of the soil wetting pattern profile at the end of infiltration are shown in Figure 6, on the basis of the following: initial soil water content θ₀ = 0.130 cm³ cm⁻³, water depth h₀ = 10.5 cm, soil bulk density ρ = 1.30 g cm⁻³, KNO₃ concentration C = 250, 600, and 900 mg L⁻¹ of clay loam.

The distribution range of NO₃⁻-N increased with the increase of KNO₃ concentration (Figure 6). At the end of infiltration (t = 360 min), the vertical distribution depth of NO₃⁻-N increased by 4.37% and 8.37% with C being 600 and 900 mg L⁻¹ compared with 250 mg L⁻¹, and the horizontal direction increased by 5.02% and 9.17%, respectively. The results indicated that the soil infiltration capacity increases with the increase of KNO₃ concentration, which is consistent with the results of Liu et al. (2017). Furthermore, it should be noted that with the increase of KNO₃ concentration, the distribution range of NO₃⁻-N gradually increases, but it does not increase linearly with KNO₃ concentration and peaks may exist, which requires further study. The results also showed that when the concentrations of KNO₃ were 250 and 600 mg L⁻¹, the concentration of NO₃⁻-N gradually increased along the direction of the wetting front and then decreased gradually; meanwhile, the concentration of NO₃⁻-N gradually decreased along the direction of the wetting front when C was 900 mg L⁻¹. The reason is that the concentration of NO₃⁻-N in the solution was 153 and 368 mg L⁻¹ when the KNO₃ concentration was 250 and 600 mg L⁻¹, both of which were lower than the initial concentration of NO₃⁻-N in the soil (496 mg L⁻¹), so that the initial NO₃⁻-N in the soil migrated to the lower layer along with water infiltration. However, when the concentration of KNO₃ was 900 mg L⁻¹, the concentration of NO₃⁻-N in the solution was 552 mg L⁻¹, which is higher than the initial concentration of NO₃⁻-N in the soil, making the concentration distribution of NO₃⁻-N similar to the soil water distribution, and is consistent with the results of Cote et al. (2003) and Ebrahimian et al. (2013).

The vertical distribution depth of K⁺ increased by only 1 cm when C was 900 mg L⁻¹ compared with 250 mg L⁻¹ (Figure 6); therefore, no obvious change exists in the distribution range of K⁺, and the water was mainly concentrated near the bottom of the furrows. The results show that the change of KNO₃ concentration has little effect on the distribution range of K⁺. Furthermore, the K⁺ concentration near the furrow bottom increased obviously with the increase of KNO₃ concentration. For example, when KNO₃ concentration was 900 mg L⁻¹, the K⁺ concentration at the bottom of the furrow increased by 248.53% compared with that at 250 mg L⁻¹. The reason may be that the K⁺ concentration increases per unit volume of infiltration with the increase of KNO₃ concentration, which further indicates that K⁺ has strong adsorption characteristics.

The distributions of NO₃⁻-N and K⁺ concentrations of soil wetting pattern profiles at the end of infiltration are shown in Figure 7, on the basis of the KNO₃ concentration C = 900 mg L⁻¹, initial soil water content θ₀ = 0.130 cm³ cm⁻³, soil bulk density ρ = 1.30 g cm⁻³, water depth h₀ = 5.5, 10.5, and 14.5 cm of clay loam.

The h₀ has little effect on the vertical distribution depth of NO₃⁻-N, but it has an obvious effect in the horizontal direction (Figure 7). At the end of infiltration (t = 360 min), when h₀ was 14.5 and 10.5 cm compared with 5.5 cm, the distribution range of NO₃⁻-N increased by 1.10% and 2.71% in the vertical direction, and in the horizontal direction increased by 20.48% and 34.87%, respectively. The reason may be that h₀ was 5.5–14.5 cm, which causes the pressure potential change to be small, so that the difference in the vertical wetting front movement distance is small. For the horizontal direction, as h₀ increased, the horizontal water surface width was increased, so that the distance between the horizontal wetting front and the center of the furrow increased, which is consistent with the results of Zhang et al. (2012). Therefore, the results suggested that h₀ should be appropriately increased to enhance the horizontal distribution of NO₃⁻-N.

The influence of h₀ on the K⁺ distribution was mainly in the horizontal direction (Figure 7). When h₀ was 14.5 cm, the horizontal distribution ranges of K⁺ were 18.49% and 26.46% higher than for 10.5 cm and 5.5 cm, respectively. The reason is that the horizontal water surface width was increased as h₀ increased, so that the distance between the horizontal wetting front and the center of the furrow increased, thereby increasing the horizontal distribution range of K⁺.

• (1)

  Soil hydraulic and solute transport parameters under different potassium nitrate (KNO₃) concentration conditions were determined by inversion solution and then verified in HYDRUS-2D. The results showed that the soil saturated water content (θ_s), empirical parameter (n), and saturated hydraulic conductivity (K_s) increase with the increase of KNO₃ concentrations, but the empirical parameter (a) shows a decrease. The simulated values with HYDRUS-2D, on the basis of soil hydraulic and solute transport parameter inversion results, closely matched the measurements of cumulative infiltration, soil water content and solute transport, and the results indicated that the soil hydraulic and solute transport parameters determined from the inversion solution were reliable.

• (2)

  The effects of initial soil water content (θ₀), potassium nitrate (KNO₃) concentration, and water depth (h₀) on the transport characteristics of nitrate nitrogen (NO₃⁻-N) and potassium (K⁺) were analyzed The results proved that the distribution range of NO₃⁻-N increases with the increase of θ₀ and KNO₃ concentration. The h₀ has little effect on the vertical distribution of NO₃⁻-N, but it has an obvious influence on the horizontal distribution, and the horizontal distribution range increases as h₀ increases. K⁺ is mainly distributed within approximately 10 cm of the bottom of the furrow, and θ₀ and the KNO₃ concentration have little effect on the K⁺ distribution range, whereas h₀ has an obvious influence on the K⁺ horizontal distribution, and the horizontal distribution range increases with the increase of h₀.

According to the results of this study, it is suggested that h₀ should be appropriately increased during furrow irrigation to increase the horizontal distribution of NO₃⁻-N. To reduce the risk of deep leaching of NO₃⁻-N, farmers should avoid fertigation when the soil water content is high in the range of suitable water contents for crop growth. Because K⁺ easily accumulates, the concentration of KNO₃ should not be too large; if it is too high, salt can accumulate on soil surfaces. Future research will require a more suitable concentration range of KNO₃ because only four concentrations were considered in this study. In addition, nitrogen mineralization, nitrification, and denitrification were neglected, which may have influenced the results of this study.

This research was supported by grants from the National Natural Science Foundation of China (No. 51579205, 51909208), Natural Science Basic Research Program of Shaanxi (No. 2019JM-063), and Scientific Research Program by the Shaanxi Provincial Education Department (20JS099).

All relevant data are included in the paper or its Supplementary Information.