Evaluating Soil Water–Salt Dynamics under Brackish Water Drip Irrigation in Greenhouses Subjected to Localized Topsoil Compaction

Abstract

Localized soil compaction in greenhouses resulting from less frequent tillage operations and frequent trampling by farmers inevitably disturbs the continuity and homogeneity of soil’s hydraulic properties, which impacts the precision of greenhouse cultivation regarding water supply and salinity control. However, predicting water–salt dynamics under partly compacted topsoil is difficult because of the interactions between many factors related to soil properties, including irrigation method and water quality, which are especially subjected to varied compaction sizes and positions. Here, two field treatments were conducted in brackish water (3 g L⁻¹) drip-irrigated plots, with the designed soil compaction region (40 cm width and 30 cm depth) adjacent to (T1) and below (T2) the drip lines. The calibrated and validated HYDRUS-2D model was applied to analyze salt exchanges across the vertical and horizontal interfaces between the compacted and non-compacted zones and the associated solute concentration variations within these two zones. The results indicated that the limited horizontal solute flux under T1 enhanced the subsequent downward flux below the drip lines, whereas, under T2, the restricted downward flux with relatively limited improved horizontal salt spreading resulted in more salt retention in the soil profile. Additional scenario simulations considering the vertical and horizontal extension of soil compaction sizes (ranging from 10 × 10 cm to 40 × 40 cm) were also conducted and revealed that, with the same increment in compaction size, the vertical extension of the compacted zone aggravated salt accumulation compared with that of horizontal extension, while the simulated cumulative water and salt downward fluxes were positive in relation to the compaction sizes in both vertical and horizontal directions under T1, but negative under T2. The findings of this study explore the effect of relative positions between drip lines and the soil compaction zone on salt transports under brackish water irrigation and reveal the potential soil salinization trend as extending compaction regions in the vertical or horizontal direction.

1. Introduction

On a global scale, intensive agriculture has been the preferred practice for shortening crop rotation periods to maximize crop yields and facilitate rapid economic benefits [1,2,3]. However, fields experiencing long-term intensive cultivation (e.g., for agronomic demand) inevitably suffer from soil compaction, which generally leads to lower structural porosity and higher structural deformation [4]. Additionally, the heavy use of farm implements enforces a high axle load and ground pressure on the surface soil and may gradually react on deeper layers [5]. It has been widely reported that soil compaction leads to soil porosity reduction with a concomitant increase in soil bulk density and strength. It also induces a decrease in hydraulic conductivity and further inhibits nutrient availability due to restricted solute convection [6]. Furthermore, the position, size, and degree of soil compaction in fields caused by human operations are often uncertain and dynamically change within the soil profile [7]; therefore, soil properties between compacted and non-compacted regions could further form capillary barriers, which ultimately reduces the water/nutrient use efficiency as a result of decreased irrigation uniformity [8]. Such negative impacts of soil compaction can be further exacerbated by the adoption of water-saving irrigation technologies, as the controlled water supply for root zones induces local-scale wetted volumes that are susceptible to variations in soil hydraulic conductivity [9].

Concerning brackish water irrigation, the salt retained by the soil profile from dripper-discharged water is the primary source for soil salinization, especially for soils with strong water-holding ability [10]. As indicated by previous research [11,12], although compaction decreases soil hydraulic conductivity at and close to saturation, it may result in improved connectivity between smaller pores and increase the conductivity for large water tensions (i.e., water contents below the field capacity). Therefore, when irrigation water contains high levels of dissolved salts, the compacted soil zone intensifies the capillary force, aggravatingly detaining the saline water flow and leading to an increased risk of salt accumulation. In addition, the water spreading pattern during the point source infiltration process and the characteristics of soil water redistribution after the cessation of water application were both found to be significantly influenced by the hydraulic properties of the soil surrounding the drippers [13]. In fields with localized soil compaction, drippers may attach to either compacted or non-compacted soil zones, thus highlighting the importance of investigating potential disparities in soil water–salt behaviors throughout an entire irrigation cycle.

The size of the compacted zone varied inconsistently along various directions within the soil profile [6], resulting in a complicated soil heterogeneity. As a result, previous studies have evaluated soil water–salt dynamics mainly on the basis of a specific or selected compacted soil situation. For example, Bohrer et al. [14] assessed the soil compaction caused by the soil respreading process during a mine-land reclamation and concluded that the soil’s resistance to water penetration varied by depth, with shallow compaction remaining unchanged but deeper compaction increasing. Besson et al. [15] investigated the salt transport process around the plough pan in lysimeters with undisturbed soil and revealed that the salt transport became more dispersive and heterogeneous due to the decreased macropore space combined with the distorted structural porosity of the compacted plough pan, which can increase the disparities and the retardation. However, a comprehensive understanding of the relationship between compaction size, water flow, and salt transport remains challenging due to the need for quantification in terms of the scope and shape of the studied compaction regions. This is primarily because downward infiltration is predominantly impeded by vertical soil stratification, while lateral water or salt spreading is expected to be more susceptible to horizontal soil property changes. Moreover, studies that quantitatively distinguish the effects of the vertical and horizontal extension of the compaction zones (i.e., the size variations in depth and width) on soil water–salt dynamics are relatively scarce.

The precise prediction of soil’s variably saturated hydraulic conductivity is a prerequisite for understanding water movement and salt transport processes in heterogeneous porous media. HYDRUS-2D model has been reportedly proven as an effective tool in evaluating strategies and technologies by numerically predicting and visibly displaying the moisture and salinity status within a modeled soil profile [16,17,18]. Taking advantage of a multitype boundary condition setup and adjustable initial condition input, HYDRUS can accurately describe the local-scale variation patterns of soil water and salinity for fields implementing micro-irrigation systems (e.g., drip, subsurface drip, sprinkler) [19,20], and can also simultaneously account for the solute concentration in irrigation water [21] or groundwater [22]. Moreover, HYDRUS provides flexible subregion selection, which is crucial and valuable in calculating the transient or cumulative flux exchanges between different subregions in the soil profile. Chen et al. [23] intentionally distinguished the root zone of a tomato/corn intercropping ecosystem using a vertical flow domain constructed in HYDRUS-2D, and analyzed the horizontal NO₃-N exchange process to reveal the competition for nutrient resources between different crops quantitatively. HYDRUS model also allows for the input of diverse sets of soil water–solute parameters, enabling the accurate consideration of soil heterogeneity resulting from field tillage practices or soil texture stratification [24,25]. In summary, the compacted zones within the studied soil profile can be assumed and divided as scale-fixed segments with specific properties and initial conditions, and the exchange flux across the edges of those segments can be recorded in a timely manner during the simulation period, which determines the feasibility of HYDRUS-2D in assessing the response of soil water–salt dynamics to localized compaction in brackish water drip-irrigated fields.

The main objectives of this research are (i) to construct two-dimensional numerical flow domains, based on HYDRUS-2D software, that can reflect drip-irrigated soil profiles with soil compaction regions below or adjacent to the drip lines; (ii) to simulate and compare the average solute concentrations between the compacted and non-compacted regions as the brackish water irrigation cycles increase; and (iii) to analyze the exchange of solute fluxes across the vertical or horizontal interfaces between compacted and non-compacted segments. Additional scenario simulations considering the size variations of compacted regions were also conducted (iv) to evaluate the response of solute concentration distribution patterns and the cumulative percolations associated with the compaction extending in both vertical and horizontal directions. It was hypothesized that altering the relative positions between soil compaction zones and drip lines would change the influence of horizontal or vertical compaction zone size variations on soil water–salt dynamics. In addition, salt accumulation would be formed along the interface between compacted and non-compacted zones because of the capillary barrier restricting water transport and retaining solute from irrigation water.

2. Materials and Methods

2.1. Field Experiments

In this study, the experimental greenhouse was located in a vegetable production region of Yangzhou city (Shatou Town, 32°17′ N, 119°29′ E) in the south-central part of Jiangsu province by the Yangtze River. This region is responsible for around half of the vegetable supply for the 2.4 million people in the urban area of Yangzhou, and more than 400 hm² of plastic-covered greenhouses are operated year-round for vegetable cultivation. The subtropical humid monsoon climate has an average annual precipitation of 1021 mm (over half of this is concentrated between May and September) and a mean annual temperature of 16.2 °C [26]. The study site featured sandy loam soil, with average proportions of sand, silt, and clay particles among the 0~200 cm soil layer, of 55.36%, 31.78%, and 12.86%, respectively. With long-term, high-frequency cultivations, a mean soil bulk density of 1.31 g cm⁻³ was formed at a depth of 0~30 cm, while the average soil bulk density of the 30~200 cm layer was 1.52 g cm⁻³. Furthermore, the groundwater (1.03~1.68 dS m⁻¹) obtained from a pumping well near the study site was the main irrigation resource.

The drip irrigation experiments considering the situation of soil compaction were conducted in pre-divided plots (2.5 m × 2.5 m = 6.25 m²) in a greenhouse area. Each plot contained an individual drip system with a typical layout pattern of “two drip lines and four crop rows” [23,27,28], in which the distance between two drip lines was 80 cm with a 35 cm dripper spacing along each drip line. Commonly, crop rows were positioned between two drip lines or coincided with one drip line. Thus, the horticulture aisle, which might typically suffer soil compaction due to farmers’ trampling, was present around the drip lines or between two drip lines. In this study, the compacted zones positioned adjacent to (T1) and below (T2) the drip lines were considered as two treatments (Figure 1). The experiment had a completely random design with three replicates for each treatment, and plots with no compacted zone but the same drip irrigation system were set as the control (CK). Before this field experiment, each designed compaction region (40 cm width, 30 cm depth) was excavated, backfilled, and compacted to a bulk density of 1.6 g cm⁻³ for both T1 and T2. All plots share the same drip irrigation schedule, which was developed based on the soil water potential in CK plots, where the values were monitored by three tensiometers positioned at 10 cm soil depth with a horizontal distance of 20 cm from the dripper in the region between two drip lines; irrigation was applied when the soil water tension was lower than −40 kPa. For each irrigation event, the brackish water was prepared by dissolving NaCl (purity at 99.8%) in freshwater (0.11~0.22 dS m⁻¹) at a designed concentration of 3 g L⁻¹, and the water was delivered through the drip line and applied for 3 h with an average dripper discharge flux of 1.55 L h⁻¹. Notably, vegetable plantings were not introduced throughout the experimental period, and the greenhouses were consistently ventilated.

2.2. Measurements and Methods

Soil samples for measuring water and salt content were extracted using soil augers (2 cm in bore diameter) before one day and after five days of each irrigation event. On the side of the two drip lines in each experimental plot, the total sixteen sampling points were located at four horizontal distances—5, 15, 25, and 35 cm away from a dripper—along the direction perpendicular to the drip line with four vertical depths (5, 15, 25, and 35 cm). Additionally, the drippers corresponding to the positions of sampling points differed for each sampling event, and the drill holes after soil sampling were all re-filled with surrounding soil to prevent the emergence of preferential flow that would affect soil water–salt movement during the next irrigation event. The gravimetric water content of the soil samples was measured through oven drying and further converted to volumetric water content based on the corresponding bulk density of the sampling position. Soil salt content was described by the electrical conductivities of a 1:5 soil/water suspension (EC_1:5) measured using a conductivity meter (DDBJ-350, Leici instrument factory of Shanghai Precision Scientific Instrument Co., Ltd., Shanghai, China). For monitoring the fluctuation of groundwater, a water level data logger (HOBO U20-001-03, Onset Computer Inc., Bourne, MA, USA) was installed at a 3 m depth in a 5 cm diameter observation well near the experimental greenhouse, and the water table depth varied within the range of 1.46~1.87 m during the study period (Figure 2).

Surface evaporation is a comparative and widely observed index that reflects the microclimate of the internal environment of a greenhouse as it is affected by light, temperature, humidity, field condition, and cultivation management [29,30,31]. However, the variations in soil surface evaporation can relatively fully reflect the greenhouse microclimate when the soil matric potential does not restrict (i.e., the soil is saturated) the soil water movement. Thus, the potential soil evaporation needs to be estimated during the study period. Before this experiment, six 2000 mL aluminum columns filled with soil samples (three with compacted soil and three with non-compacted soil) were gradually saturated from bottom to top using Mariotte bottles. Then, all columns buried in the soil to the top were flush with the field surface. After that, the water level reduction in each Mariotte bottle was recorded every 3 h to calculate the water consumption amount (as an approximation of potential evaporation) until the point when the Mariotte bottle could no longer supply water. In the meantime, a standard 20 cm diameter evaporation pan (DF-AM3, Oriental Xinhong Science and Technology Ltd., Beijing, China) filled with distilled water was placed at the center of the studied greenhouse to test the pan evaporation. The relationship between the pan evaporation and the amount of water consumption of the soil samples was developed (Figure 3) and subsequently utilized to estimate the soil potential evaporation during the experimental period based on the pan evaporation data (Figure 2).

2.3. Numerical Model Description and Construction

2.3.1. Governing Equations for Water Flow and Solute Transport

Soil water flow in the experimental plots under different soil compaction treatments was modeled using HYDRUS model with the software version of 5.03, which is based on the mass conservative formulation [32] with the application of a Galerkin finite-element method to solve Richards’ equation [33] for variably saturated flow in porous media. The two-dimensional Richards’ equation governing the water flow in a homogeneous and isotropic soil can be defined as follows:

$$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial x}\left(K(h)\frac{\partial h}{\partial x}\right)+\frac{\partial h}{\partial z}\left(K\left(h\right)\left(\frac{\partial h}{\partial z}+1\right)\right)-S(h)$$

(1)

where x and z are the horizontal and vertical coordinates (cm), respectively. θ is the soil volumetric water content (cm³ cm⁻³), h is the soil water pressure head (cm), K is the hydraulic conductivity (cm h⁻¹), and S is a distributed sink term related to root water uptake of plants (h⁻¹) (which was ignored in this model study).

The unsaturated soil hydraulic properties described rely on the application of a soil water retention curve and soil hydraulic conductivity function in the formulation of a van Genuchten–Mualem model [34,35]:

$$\theta (h)=\left\{\begin{array}{c}{\theta }_r+\frac{{\theta }_s-{\theta }_r}{{\left[1+{\left(\alpha h\right)}^n\right]}^m},\left(h<0\right)\\ {\theta }_s,\left(h\ge 0\right)\end{array}\right.$$

(2)

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

(3)

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

(4)

where θ_s and θ_r are the saturated and residual water content (cm³ cm⁻³), respectively, K_s is the saturated hydraulic conductivity (cm h⁻¹), S_e is the relative saturation (dimensionless), α (cm⁻¹), n, and m (dimensionless) are shape parameters with m = 1 − 1/n, and l is a pore connectivity parameter with an estimated value of 0.5 for average conditions in a range of soils [34].

The advection–dispersion equation adopted by HYDRUS-2D was computed for describing the solute transport in a homogeneous soil, neglecting chemical reactions and interactions with the solid phase. The two-dimensional form of this equation is written as follows:

$$\frac{\partial (\theta C)}{\partial t}=\frac{\partial }{\partial x}\left(\theta D_{xx}\frac{\partial C}{\partial x}+\theta D_{xz}\frac{\partial C}{\partial z}-q_{x}C\right)+\frac{\partial }{\partial z}\left(\theta D_{zz}\frac{\partial C}{\partial x}+\theta D_{zx}\frac{\partial C}{\partial z}-q_{z}C\right)$$

(5)

where C is the total solute concentration (mg cm⁻³), q_x and q_z are the horizontal and vertical volumetric fluid fluxes (cm h⁻¹), respectively, and D_xx, D_xz, D_zz, and D_zx are the components of the dispersion coefficient (cm² h⁻¹).

2.3.2. Modeling Domain Construction and Boundary Condition Setup

The constructed two-dimensional flow domain corresponds to the soil profile (200 cm deep and 200 cm wide) directly below the two dippers and perpendicular to the two drip lines of the experimental plot (Figure 4). The flow domain follows an axisymmetric design with the drippers of the left and right drip lines located horizontally 60 cm away from the left and right sides of the domain boundary, respectively. More than 60,000 unstructured triangular-finite-element meshes spatially discretized the flow domain, with finer mesh adopted close to the soil surface where fluxes of irrigation or evaporation occurred. Furthermore, the 0~30 cm soil layer of the flow domain was uniformly divided into five (30 cm deep and 40 cm wide) rectangle subregions (Figure 4) for the representation of localized soil compaction; i.e., T1 specifies the compacted soil in subregions 1, 3, and 5 and non-compacted soil in subregions 2 and 4, while T2 corresponds to the opposite pattern. The remaining region of the domain (30~200 cm) was defined as subregion 6 and set up with the same soil condition for all treatments. Meanwhile, observation points for recording the dynamic variability in the modeled soil water–salt content were set in this flow domain and followed the sixteen sampling positions of the field measurements.

The potential soil evaporation estimated by pan evaporation (Figure 2 and Figure 3) was imposed as the atmosphere condition boundary for the surface of this flow domain, and the potential transpiration was ignored as no plants were involved in this experiment. Meanwhile, the lateral flow between plots (average spacing of 1.5 m) was ignored and the no-flux boundary condition was specified on the right and left sides of the flow domain. At the bottom of the domain geometry, a variable pressure head with a constant concentration of 1.42 mg cm⁻³ was imposed for maintaining the measured water table depths and average groundwater salinity. Furthermore, the line-source infiltration assumption [36,37] was adopted to represent the intermittent flow from the dripper, and two variable-flux boundaries were assigned below the locations of the two drippers (Figure 4). An approximate relationship expressing a disc-source radius (R_s), as proposed by Wooding [38] and Bresler [39], was applied to estimate the width of the variable-flux line:

$$R_s={\left(\frac{4}{{\alpha }^2{\pi }^2}+\frac{q}{\pi K_s}\right)}^{1/2}-\frac{2}{\alpha \pi }$$

(6)

where K_s and α refer to the soil hydraulic parameters in Equations (2) and (3), and q is the dripper discharge rate (cm³ h⁻¹); this is assumed to be a steady-state flow during each irrigation event. Since the vertical model domain is at the dripper position and perpendicular to the direction of the drip lines, the length of a variable-flux boundary was set as 2R_s on the upper boundary for each dripper.

2.3.3. Initial Conditions and Parameter Optimization

The soil water content and salt concentration measured one day prior to the first irrigation event were assumed as the initial conditions. Soil samples of the initial condition measurements, accounting for the compacted and non-compacted soil regions, were collected every 15 cm down to a depth of 180 cm, and the initial distributions of soil water content and salt concentration were assumed to be uniform by setting the corresponding average measured values within each subregion. Additionally, three kinds of soils with different bulk densities were considered in this model study, i.e., non-compacted soil from 0~30 cm depth with 1.31 mg cm⁻³ bulk density; compacted soil from 0~30 cm depth with 1.60 mg cm⁻³ bulk density; and deep soil from 30~200 cm depth with 1.52 mg cm⁻³ bulk density. During the model operation, the salt content in the liquid phase was considered for the quantitative analysis of salt mass flux during the simulation. Therefore, the measured values of soil EC_1:5 (dS m⁻¹) were converted to the salinity in the actual soil water content according to the following relationships:

$$E{C}_e=8.02\times E{C}_{1:5}+0.23$$

(7)

$$E{C}_{sw}=\frac{E{C}_e\times SP}{\theta }$$

(8)

$$TDS=0.681\times E{C}_{sw}$$

(9)

where EC_e is the electrical conductivity of the saturated paste extract (dS m⁻¹), EC_sw is the electrical conductivity of the actual soil water content of θ, SP is the saturation percentage %, and TDS is the total dissolved solids (mg cm⁻³) in the soil extract and is used to represent the solute concentration in the liquid phase during the model study. The linear relationships of Equations (7) and (9) were developed in laboratory, and Equation (8) was established according to another HYDRUS model study examining water and salt dynamics [40].

Soil hydraulic parameters, as described by the van Genuchten–Mualem model (i.e., θ_r, θ_s, α, n, and K_s in Equations (3) and (4)) were initially estimated using the pedotransfer function based on the neural networks [41], which is based on the input information of soil particle size and bulk density. Previous studies reported that n, θ_s, and K_s are the three most sensitive parameters for soil water flow [27]. Thus, the laboratory-tested values of θ_s and Ks before the experiment and ROSETTA-predicted n were set as the initial values and subsequently modified for parameter optimization, while the less sensitive parameters of θ_r and α were fixed at their initial values for all the three soils during the parameters’ optimization. In terms of solute transport parameters, soil longitudinal (D_L) and transverse dispersivities (D_T) were known as scale-dependent parameters [42], which variated by the spatial scale of field or laboratory measurements [16]. In this study, D_L was initially set as 20 cm, which is suggested to be one-tenth of the depth of the flow domain [43,44,45], and the ratio of D_L/D_T = 8 was selected from the literature related to the texture of sandy loam [46]. The molecular diffusion coefficient was fixed at 1.63 cm² day⁻¹ following a previous model study in which the author used NaCl as the primary solute [24].

Parameter optimization was conducted with calibration and validation procedures. During the calibration, the soil hydraulic and solute parameters of non-compacted and deep-layer (30~200 cm depth) soils were optimized based on the comparison between the measured and modeled data of soil water–salt content under CK using the trial–error approach until the simulation accuracy requirements were met. Further, the optimized parameters of non-compacted and deep-layer soils were fixed and implemented into the model domain of T1 for calibrating the corresponding parameters of compacted soil. The validation was performed on the basis of the measurements of T2 to verify the applicability and accuracy of the calibrated soil parameters.

2.4. Statistical Analysis

The agreement between the field-measured and model-simulated data for reflecting the simulation accuracy of parameter calibration and validation was evaluated using the statistical indicators of the root-mean-square error (RMSE), the Nash–Sutcliffe model efficiency coefficient (NSE), and the mean relative error (MRE):

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{\left(M_i-S_i\right)}^2}}$$

(10)

$$NSE=1-\frac{{\displaystyle {\sum }_{i=1}^n{\left(M_i-S_i\right)}^2}}{{\displaystyle {\sum }_{i=1}^n{\left(M_i-\overline{M}\right)}^2}}$$

(11)

$$MRE=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n\frac{\left|S_i-M_i\right|}{S_i}}\times 100\%$$

(12)

where M and S refer to the measured and simulated data of soil volumetric water content or EC_1:5, respectively. $\overline{M}$ is the mean value of the measured data, and n is the number of observations. RMSE and MRE values closer to zero indicate a better simulation performance. The NSE generally varied between −∞ and 1.0, with the best fitting being NSE equal to 1; an acceptable level of model performance was attained when the NSE values were in the range of 0 to 1 [43].

2.5. Modeling Scenarios

Soil regions with increased bulk density caused by mechanical loading are likely to extend compaction cycles. Therefore, numerical simulations based on the parameter-optimized 2D model were further adjusted and operated to explore the effect of compaction size on soil salt dynamics. In Figure 5, sixteen scenarios were designed, taking into account the size extension of the compaction region along vertical (10, 20, 30, and 40 cm depths) and horizontal (10, 20, 30, and 40 cm widths) directions from the position below the drippers, and were implemented for both T1 and T2 conditions. Among the designed scenarios, the S43 scenario was equal to the field experimental conditions and is also included for comparison with the other scenarios.

3. Results

3.1. Model Calibration and Validation

The field-measured volumetric soil water contents at different sampling locations during the experimental period were graphically compared with the corresponding simulated values of the observation points (Figure 6). Under three treatments, the simulated values generally displayed good agreement with the corresponding measured values, in that the variation ranges of soil water content were enlarged with increasing lateral distance from the dripper, and relatively higher water content was obtained in deeper layers. Meanwhile, the established 2D model more fairly approximated the measured data related to the compacted region than those for non-compacted regions, as shown by the relatively lower RMSE and MRE and higher NSE (

Table 1. Statistical comparison of soil water content (SWC) and EC_1:5 at 5, 15, 25, and 35 cm lateral distance from the dripper within 0~35 cm soil depth during calibration (CK, T1) and validation (T2).

Treatment	Lateral Distances from the Dripper	Soil Water Content			EC1:5
		RMSE	NSE	MRE	RMSE	NSE	MRE
CK	5 cm	0.018	0.767	8.202	0.057	0.650	7.245
	15 cm	0.023	0.714	7.416	0.067	0.667	5.531
	25 cm	0.019	0.732	7.510	0.091	0.671	5.969
	35 cm	0.017	0.743	6.989	0.132	0.577	9.873
T1	5 cm	0.019	0.721	6.377	0.066	0.641	7.051
	15 cm	0.024	0.805	5.417	0.065	0.710	7.120
	25 cm	0.016	0.654	6.072	0.087	0.539	8.781
	35 cm	0.017	0.719	6.085	0.090	0.531	8.159
T2	5 cm	0.024	0.729	5.349	0.053	0.559	6.217
	15 cm	0.017	0.634	5.908	0.064	0.602	6.916
	25 cm	0.026	0.652	6.277	0.078	0.562	6.562
	35 cm	0.025	0.673	5.883	0.101	0.542	8.249

Note: Data are the means of five replications, and those followed without the same letter differ significantly at the p = 0.05 level. Soil textures were determined using the USDA (United States Department of Agriculture) textural soil classification system.

). The RMSE, NSE, and MRE values at different lateral distances varied from 0.017 to 0.026 cm³ cm⁻³, from 0.652 to 0.805, and from 5.417 to 8.202%, respectively, for the non-compacted region with the three treatments. For the compacted regions of T1 and T2, these values were from 0.016 to 0.022 cm³ cm^-3, from 0.634 to 0.729, and from 5.349 to 6.085 %, respectively. Regarding the soil salt simulation, the model-simulated values of liquid-phase solute concentration were converted to EC_1:5 (based on Equations (6) and (7)) for comparison with the measured values (Figure 7). The dynamics of the measured and simulated salt contents of the three treatments followed similar variation patterns; the values were likely to be varied within a high and wide range within a shallow depth and the increasing trend in salt content was possibly more pronounced with the increase in lateral distances from the dripper. The statistical comparison in the salt simulation showed a wider error range than the soil water simulation (Table 1). The RMSE, NSE, and MAE values obtained by comparing measured and simulated soil EC_1:5 for three treatments at different lateral distances to the dripper were in the ranges of 0.053 to 0.132 cm³ cm⁻³, 0.531 to 0.710, and 5.531 to 9.873%, respectively. However, lower NSE values and higher RMSE and MAE values were recorded with increased lateral distance from the dripper, while the statistical indicators presented no intuitive increase or decrease between the compacted and non-compacted regions. Overall, the statistical errors between field-measured and HYDRUS-simulated values were within an acceptable limit during the calibration and validation phases, suggesting that the established model is capable of accurately predicting the water–salt dynamics in drip-irrigated fields with soil compaction. The initial and verified model input parameters are listed in

Table 2. Soil parameters for water flow and solute transport in the established HYDRUS-2D model.

	Bulk Density (g·cm−3)	Soil Hydraulic Parameters						Solute Transport Parameters
		θr (cm3·cm−3)	θs (cm3·cm−3)	a (cm−1)	n	Ks (cm·h−1)	l	DL (cm)	DT (cm)
Initial parameters	1.31	0.049	0.452	0.016	1.475	5.25	0.5	20	2.5
	1.52	0.045	0.403	0.021	1.427	3.43	0.5	20	2.5
	1.6	0.042	0.396	0.025	1.378	2.27	0.5	20	2.5
Optimized parameters	1.31	0.049	0.44	0.016	1.48	3.50	0.5	60	7.5
	1.52	0.045	0.415	0.021	1.50	2.20	0.5	45	5.63
	1.6	0.042	0.36	0.025	1.32	1.75	0.5	35	4.38

.

3.2. Soil Salt Distribution and Dynamic Transport Pattern

The axisymmetric modeling domain design results in the same values of mean salt concentration (MSC) in subregions 1 and 5 and in subregions 2 and 4 for all treatments (Figure 8). In subregions 2 and 4, the MSC below the drip line was expected to be greatly decreased after the irrigation leaching, and was further slowly raised after eight irrigation applications of brackish water. In contrast, the MSC values of subregion 1 or 5 under the three treatments (16.14~18.7 mg cm⁻³) gradually increased and remained almost twice as high as the corresponding values of subregion 2 or 4 (7.68~10.5 mg cm⁻³) after drip irrigation. The final MSC values (at DAS 120) in subregion 3 showed a relatively stable variation, with an 8.93% decrease under T1 and 5.82% and 10.17% increases under CK and T2 compared with the initial values, respectively. Furthermore, the final MSC values for CK in subregion 1 or 5 and subregion 2 or 4 were 1.34% and 7.02% higher than T1, but 1.92% and 1.14% lower than T2, suggesting that T1 presented the lowest salt buildup under the same brackish water irrigation strategy.

The soil salt exchanges in the vertical and horizontal directions between the compacted and non-compacted regions were captured on the interface between different subregions of the modeling flow domain (Figure 9). The exchange flux at the different interfaces followed a similar variation pattern during each irrigation event, and the process during the sixth irrigation event was selected to graphically reveal the discrepancies between treatments. The salt flux at two vertical interfaces (H1 and H2) presents a couple of symmetric variation patterns; the salt was immediately transported from subregions 1 and 3 to subregion 2 after the irrigation started, and then the direction of salt flux was gradually reversed and it was transported from subregion 2 to subregions 1 and 3. Finally, the flux approached zero as the irrigation stopped. For each irrigation event, the most pronounced flux exchanges were obtained in CK (−87.27~39.44 mg cm⁻¹ h⁻¹ at H1), while the lowest exchange of salt flux occurred in T1 (−57.43~25.64 mg cm⁻¹ h⁻¹ at H1). Regarding the salt exchange through the horizontal interfaces (V1, V2, and V3), there was no flux direction reversion like that which appeared in the vertical interfaces, and the magnitude of downward salt fluxes far exceeded the upward salt flux during the no-irrigation period. There was a considerable increase in the downward salt flux at the V2 interface, especially for T1, where the associated downward salt fluxes were 22.45 and 13.52 times larger than those at V1 and V3, respectively. Additionally, in CK and T2, the same compaction conditions of soil above and below the interface of V1 and V3 resulted in promoting salt exchange, with the average peak fluxes of CK and T2 being 78.21% and 49.36% at V1, and 75.83% and 44.03% at V2 compared with those of T1, respectively.

3.3. Soil Salt Distributions and Transport under Different Compaction Scenarios

Figure 10 shows the simulated salt distributions under T1 and T2 with the scenarios of minimum (S11) and maximum (S44) soil compaction size after the last irrigation event (DAS 108 to 120). The small compaction region of S11 (10 × 10 cm) resulted in almost the same distribution characteristics between T1 and T2 (Figure 10a,b), in that the salt concentrations were generally low after irrigation leaching (DAS 108) and the evaporation-induced increases in salt concentration were progressively more pronounced with salt concentration gradients away from the dripper. In comparison, under S44 (Figure 10c,d), there was a distinguishable relatively higher concentration region in T2 formed in the middle of the 2D flow domain compared with T1, demonstrating that the larger size of compacted soil below the drippers limited the leaching effects and induced salt accumulation downward to the soil profile.

To quantitatively characterize the salt concentration gradients caused by drip irrigation with brackish water under different soil compaction scenarios, two observation regions within a 0 to 40 cm soil depth representing the regions outside and between the two drip lines were delineated from the 2D domain; the proportion of areas with different concentration–gradient ranges from the entire observation region was calculated based on the solute concentration distribution at the end of the simulation. For regions outside the drip lines (Figure 11), the solute concentrations of different scenarios were mainly distributed in the range of 10~20 mg cm⁻³, and the proportion of higher solute concentration areas was positively correlated with the size of the compaction region. However, the proportions of high salt distributions showed a more pronounced increase as the compacted range extended along the horizontal direction than along the vertical direction for T1; however, for T2 the trend was reversed. For example, the proportions of 20~30 mg cm⁻³ solute concentration areas were increased by 3.7% and 0.9% with the horizontal extension of the compacted region from S41 to S44 in T1 and T2, respectively, while the corresponding increases with the compacted region extending vertically from S14 to S44 were 2.9% for T1 and 2.7% for T2. In addition, the proportions of 0~10 mg cm⁻³ salt distributions varied in a wider range for T2 (23.3~30.4%) than for T1 (28.9~30.8%) under different scenarios. Comparatively, in the region between two drip lines (Figure 12), there was no proportion of the area with solute concentrations over 20 mg cm⁻³, and around 70% of this region exhibited a 0~10 mg cm⁻³ solute concentration under all the simulated scenarios for both T1 and T2, suggesting the wetted volumes caused by drip irrigations were overlapped in the region between two drip lines that enhanced the water supply and associated salt leaching. Furthermore, areas with solute concentration above 15 mg cm⁻³ were obtained in the scenarios with 30 and 40 cm compaction depths under T1, whereas such areas only appeared in the scenarios with a 40 cm compacted depth under T2. Although the higher solute concentration areas were more likely to appear in T1, the proportions of the area of all the scenarios with solute concentration below 10 mg cm⁻³ under T1 (averaged at 73.83%) were mostly higher than under T2 (averaged at 71.15%).

4. Discussion

4.1. HYDRUS-2D Performance

Based on HYDRUS, two-dimensional numerical simulations of soil water–salt dynamics for drip-irrigated soil profiles are widely accepted and operated under various soil physical properties, agronomic operations, and crop rotation conditions [37,47,48]. In this study, the parametrized (soil hydraulic and solute transport) model showed a reasonable performance in representing the distribution patterns and trends of field-measured water–salt content during drip irrigation cycles, and the statistical indicators illustrated that fewer modeling errors (i.e., lower RMSE values and higher NSE and MRE values) were obtained in the simulation of soil water content than salt content during the calibration and validation procedures. A similar phenomenon was also reported by previous HYDRUS-based studies relevant to the simulation of salt dynamics under one- [49], two- [16], and three-dimensional [50] modeling domains. The possible explanation for this is that HYDRUS establishes the convection dispersion equation based on the solution of Richards’ equation and, thus, the errors from soil water simulation are reflected in the salt simulation. Additionally, frequent wetting–drying cycles during drip irrigation events make soil susceptible to the shrinking–swelling process [51], which promotes the formation of desiccation cracks (present in some experimental plots) along the interface between compacted and non-compacted zones, and inevitably induces some random and localized changes in hydraulic porosity. Thus, setting fixed hydraulic and solute parameters seems to be a potential source of statistical errors in modeling soil water–salt content for a drip-irrigated soil profile, especially when soils with different bulk densities are alternately distributed.

4.2. Salt Dynamics among Compacted and Non-Compacted Regions

The pulse and low flow discharge of drip irrigation make the extension of wetted areas sensitive to soil hydraulic porosities [52]. There are inherent problems regarding salt accumulation due to limited soil leaching and persistent evaporation, which aggravate soil salinity buildup, particularly at the surface and periphery of the point source-induced wetted bulbs [17,27]. The above-mentioned features have been approximately verified by these modeling results, as the soil regions below (subregions 2 and 4) and outside (subregions 1 and 5) the two drippers reflect the lowest and highest solute concentrations, respectively (Figure 8), while the solute concentration in the region between the two drippers (subregion 3) only presents a slight increase due to the enhanced soil leaching caused by wetted volumes overlapping [53]. Berezniak et al. [47] modeled drip irrigation under manipulated soil texture conditions and indicated that although the limited discharge of drippers weakens the water spreading in the soil profile, the vertical preferential flow pathway caused by the capillary barriers among soil profiles can still lead to a risk of salinization for deeper layers. Similar results were obtained in the salt flux simulation (Figure 9), in which the highest downward solute flux was obtained in the non-compacted zone of T1 due to the fact that the lateral water/salt diffusion was restricted by the vertical capillary barrier formed by the left- and right-side compaction zones. Notably, the local salt accumulation by irrigated brackish water leads to high concentration gradients along the horizontal direction and a dispersive solute flux component opposite to the direction of the water flux (i.e., H1 and H2 in Figure 9). Such phenomena are physically wrong because it appeared that molecules or ions were moving against the flow direction through dispersion; this can be ascribed to the basic theory of the convection–dispersion equation in that macroscopic dispersion is described as an isotropic random transport process [54]. By contrast, such phenomena were not significant along the horizontal interface (i.e., V1, V2, and V3 in Figure 9) because the gravitational potential strengthened the downward water flow and the concentration gradient was relatively lower along the vertical direction.

Moreover, in line with some previous field or model studies on drip irrigation with saline water [48,55], shallow groundwater (varying within a 1~3 m depth) probably plays an insignificant role in replenishing the salt rise caused by evaporation. As observed in this study (Figure 9), the downward solute flux induced by irrigation from the upper layers far exceeds the upward flux induced by evaporation, indicating that the salt accumulation in the soil profile primarily relies on the storage of salt from irrigation water.

4.3. Effect of Compaction Sizes on Salt Distribution

The capillary barrier formed due to hydraulic property discrepancy may significantly prevent the solute’s return or release across different interfaces, primarily when the soil matrix dominates the capillary-driven advection fluxes after irrigation [47]. The results of the scenario simulations indicated that the solute concentration, within 0~40 cm of soil depth, generally increased with the size of the compacted region (Figure 11 and Figure 12), while the vertical size extension of the compaction region led to a higher rise in solute concentration compared with an equivalent extension in the horizontal direction. The possible explanation for this is that the vertical capillary barrier obstructs the horizontal solute flux induced by the combined effect of water pressure and capillary suction during drip irrigation, while the vertical flux is also facilitated by the gravitational potential, which partially counteracts the obstruction effect of the horizontal capillary barrier. Furthermore, the semi-elliptic cylindrical wet bulb formed by dripper discharge reportedly expands more widely but less deeply for soil with a lower hydraulic conductivity [56]. This shape response of the wetting bulb can further elucidate why the impact of vertical capillary barrier size on solute concentration is more pronounced under T2 conditions, where lateral water–salt movements can be further intensified due to compacted regions being located directly below the drippers. However, this is inconsistent with previous two-dimensional simulation studies [24,57] that observed noticeable disruptions in salt distribution on the interface between soil layers with different hydraulic properties. There was no sharp increase or decrease in solute concentration at the periphery of the compacted region regardless of the size of the compaction region (Figure 10), which was probably due to the uniform soil texture in the soil profile that offset the impact of the capillary barrier caused by the difference in bulk density on salt migration and accumulation. For both T1 and T2, the soil region between two drip lines that received irrigation water from two drippers showed a distinctly lower salinity level than the region outside the drip lines, implying that the intensity of soil leaching plays a crucial role in creating low-salinity areas in a soil profile when adopting brackish water for drip irrigation.

It should be noted that, in establishing a model domain for different scenarios, we hypothesized that the soil was evenly compacted into a high-bulk-density rectangular region, resulting in homogeneous soil hydraulic properties for water and salt movement in both compacted and non-compacted regions. However, in realistic conditions, the size and shape of soil compaction zones are difficult to accurately predict because of the interactions between many factors related to soil properties, cultivation practices, and field management, which lead to irregular alternations in pore size distribution and associated fluctuations in bulk density [6]. Therefore, the simplifications in the field experiment and numerical simulations of this study only allow for an approximate assessment of the influence of the horizontal or vertical extension of the compaction zone on soil water and salt transports, limiting the accuracy of the results. Therefore, it is necessary to further demonstrate the relevant conclusions in light of situations where there are significant nonlinear variations or random distributions in bulk density within the compaction zone.

5. Conclusions

The present study tried to understand the impact of localized topsoil compaction on salt transport and distribution under drip irrigation with brackish water. Two case studies corresponding to soil profiles with compacted segments below and adjacent to the drip lines were experimentally and numerically (HYDRUS-2D) studied to assess the response of soil salt dynamics to the relative positions between drip lines and soil compaction zones. The variation in simulated solute flux across the vertical and horizontal interfaces between compacted and non-compacted regions suggested that a pronounced capillary barrier effect of the vertical interface was obtained as the salt migrated from the non-compacted region to the compacted region. Meanwhile, the irrigation-induced downward solute fluxes below the drippers were further enhanced as the compacted region was located at the two sides of the drip lines. In contrast, the magnitude of upward solute flux during the no-irrigation period could almost be neglected compared with the downward solute fluxes for all the compacted conditions. The results of the exploratory scenario simulation indicate that the vertical size extension of the compaction region led to a higher rise in soil salt content compared with an equivalent size extension along the horizontal direction, which means that limiting the vertical development of the compaction zone helps to mitigate the salt accumulation under brackish water drip irrigation. Furthermore, the enhanced soil leaching induced by the overlapped wetting volumes between the two drip lines could offset the soil salinization aggravated by the extension of soil compaction. In terms of practical applications for greenhouse farmers to improve their soil management practices, it is suggested that (1) horticulture aisles for farmers’ walking should be widened to avoid the vertical extension of the compacted zone caused by high-frequency trampling over a narrow area; and (2) crops should preferably be planted between two drip lines when the wetted volumes of two drip lines could overlap.