Development and validation of a model for soil wetting geometry under Moistube Irrigation

We developed an empirical soil wetting geometry model for silty clay loam and coarse sand soils under a semi-permeable porous wall line source Moistube Irrigation (MTI) lateral irrigation. The model was developed to simulate vertical and lateral soil water movement using the Buckingham pi (π) theorem. This study was premised on a hypothesis that soil hydraulic properties influence soil water movement under MTI. Two independent, but similar experiments, were conducted to calibrate and validate the model using MTI lateral placed at a depth of 0.2 m below the soil surface in a soil bin with a continuous water supply (150 kPa). Soil water content was measured every 5 min for 100 h using MPS-2 sensors. Model calibration showed that soil texture influenced water movement (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p$$\end{document}p < 0.05) and showed a good fit for wetted widths and depths for both soils (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$nRMSE$$\end{document}nRMSE = 0.5–10%; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$NSE \ge$$\end{document}NSE≥ 0.50; and d-index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ge$$\end{document}≥ 0.50. The percentage bias \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( {PBIAS} \right)$$\end{document}PBIAS statistic revealed that the models’ under-estimated wetted depth after 24 h by 21.9% and 3.9% for silty clay loam and sandy soil, respectively. Sensitivity analysis revealed agreeable models’ performance values. This implies the model's applicability for estimating wetted distances for an MTI lateral placed at 0.2 m and MTI operating pressure of 150 kPa. We concluded that the models are prescriptive and should be used to estimate wetting geometries for conditions under which they were developed. Further experimentation under varying scenarios for which MTI would be used, including field conditions, is needed to further validate the model and establish robustness. MTI wetting geometry informs placement depth for optimal irrigation water usage.

Agriculture is the largest consumer of blue water 1 at 70% of all global freshwater resources 2,3 . Novel irrigation technologies such as sub-surface irrigation and porous pipes promote water conservation 2 . Moistube Irrigation (MTI) is a semi-permeable porous pipe that has reported improved field water use efficiency (fWUE). MTI is a sub-surface irrigation technology whose discharge is facilitated by an applied pressure, or at zero pressure, it utilizes soil water matric potential ( ψ ) that causes a pull effect, thus facilitating discharge. MTI is a semi-permeable porous wall tubing; thus, its wetting geometry is classified as a line source emitter 4,5 . Soil water movement under various irrigation technologies has informed irrigators on the effective placement depth and lateral spacing that promote crop water use efficiency (cWUE) and fWUE.
There is limited empirical knowledge on models that facilitate the estimation of wetted perimeters under porous wall emitters. Knowledge of soil wetting geometry is critical in optimizing MTI irrigation system design (lateral placement depth and spacing) and operation (discharge rates, irrigation set times and satisfying irrigation water requirements). To maximize the advantages offered by sub-surface irrigation, knowledge of soil wetting geometries aids in irrigation network design, i.e., emitter spacing and placement depths, which subsequently improve irrigation schedules 6 , minimize run-off losses, promotes higher irrigation uniformity 7,8 , increases water productivity (WP) and fWUE 4,9 . Soil wetting geometries can be determined either experimentally or using modelling tools. The former is expensive and time-consuming. Modelling is a time-saving exercise, and numerical models have gained wide applicability over their counterparts (analytical and empirical models) because of their robustness and use of Using t Buckingham π theorem 16,19 , four πs were derived as presented in Eqs. (3)- (7). The four πs were derived because the Buckingham π theorem states that if there is a physically meaningful equation involving a certain (1) (W, Z) = f (V , q, k, D) (2) f (V , q, k, W, Z, D) www.nature.com/scientificreports/ number, n, of physical variables in a problem and these variables contain 'm' primary dimensions, the equation relating all variables will have (n-m) dimensionless groups. There were six variables with two primary variables, namely W and Z , which resulted in 4 dimensionless groups.
where Multiplying the π 3 and π 4 yielded the dimensionless soil water content per unit length of MTI (V * ) as presented by Eq. (8).
By taking the square root of the product of π 4 and (π 2 ) 2 yielded the dimensionless wetted width (W * ) as presented by Eq. (9).
The square root of the product of π 4 and (π 1 ) 2 yielded the dimensionless wetted depth (Z * ) as presented by Eq. (10).
Schwartzman and Zur 15 and Singh, Rajput 19 postulated that there exists a relationship amongst dimensionless parameters. For this research, the relationships are as presented in Eqs. (11) and (12); where A 1 , A 2 , b 1 , and b 2 are constants for a 2-dimensional flow model. The constants A 1 and b 1 , were determined from the graphical plot of V * and W * whereas the constants A 2 and b 2 were determined from the graphical plot of V * and Z * .
Combining Eqs. (8) and (9) and Eqs. (8) and (11) yielded the wetted width ( W ) and wetted depth ( Z ) functions presented in Eqs. (13) and (14), respectively. Experimental design and data collection. Soil hydraulic parameters and textural characteristics. The silty clay loam (34% clay, 58% silt, 8% sand) was obtained from the University of KwaZulu-Natal's Ukulinga Research Farm in Pietermaritzburg, KwaZulu-Natal, South Africa (29° 39′ 44.8ʺ S 30° 24′ 18.2ʺ E, altitude: 636 m). The coarse sand soil (98% sand and 2% gravel) was obtained from Genie sand in Pinetown, KwaZulu-Natal, South Africa (29° 48′ 08.7ʺ S 31° 00′ 37.8ʺ E). Soil samples were subjected to soil textural analyses using the hydrometer method. The experiment sampled five depths for textural analysis, and the resultant textural data was fed into the SPAW model (Saxton and Willey, 2005) to determine saturated hydraulic conductivity ( k s ) was derived. Other soil hydraulic parameters total porosity ( θ r ), residual soil water content ( θ s ), and shape fitting parameters ( n , m , and α) ( Table 1) were laboratory determined using the soil-water retention pressure method 4,29,30 . The 50 cm depth soil sample for the silty clay loam was used to fit the van Genuchten parameters (3) f (π 1 , π 2 , π 3 , π 4 ) = 0 www.nature.com/scientificreports/ because the 50 cm plot provided a smooth curvilinear shape and the resultant parameters closely aligned with Rawls, Brakensiek 31 . The sandy soil was commercially acquired hence the absence of varied sampling depths. The methods were selected based on the reliability of results and equipment availability.
Measurement of soil wetted front. The soil was air-dried, crushed and sieved through a 2 mm sieve. Thereafter, the soil was loaded into a soil bin measuring 1 m (H) × 1 m (W) × 0.5 m (B). The soil bin had transparent Plexiglass walls, and soil loading was done gently to avoid compaction and possible crushing of the MTI tubing. To prevent MTI collapse under the soil surcharge, MTI was supplied with water before loading the soil till it was turgid. To prevent MTI smearing and potential nano-pore blocking, the water was supplied upon reaching the MTI burying level. The MTI lateral was placed at a depth of 0.2 m below the soil surface, and upon soil loading, MPS-2 sensors were simultaneously installed at prescribed depths ( Table 2). The initial soil-water content for the silty clay loam was 1.02 × 10 -6 m 3 and 6.64 × 10 -7 m 3 for the sandy soil. Both soils were packed at a bulk density of 1.4 g cm −3 . The MPS-2 sensors measured water potential (− 10 to − 500 kPa) and temperature, and they were calibrated by soaking them in de-ionised water for a period of 72 h before installation. The de-ionised water was used to substitute for the conventional mercury porosimeter experiment for calibration. The calibration curve was obtained from the MPS-2 and MPS-6 operators manual (Decagon Devices Inc, 2017). Water to the MTI lateral was supplied at a pressure head of 150 kPa, which discharged 2.39 l h −1 m −1 ( Table 3).
The experiment was carried out in two phases. The first dataset of measured variables was used for model calibration, whilst the measured variables from the second phase dataset were used for model validation. The measured variables are summarised in Table 4. Both the first and second phases were carried under identical conditions. Soil water-retention curves derived from the soil-water retention experiment were used to determine the volumetric soil water content. The experimental equipment set-up is shown in Fig. 1.  Table 2. MPS-2 sensors placement depths and lateral spacing for the respective soils.   process was done to determine the models' constants A 1 , A 2 , n 1 and n 2 . To improve the models' precision and accuracy, iterations were carried out on constants A 1 and A 2 .

Models' validation.
The parameter variability-sensitivity analysis (PV-SA) method was employed to validate the models (Sargent, 2013). A sensitivity analysis was carried out to assess the effects of k , q , and D on both W and Z . The sensitivity analysis was done by holding all the other variables constant and assessing the functional relationship a particular "active variable" had with W and Z . Table 4 summarise the two independent experimental phases and their respective purposes.

Models' evaluation.
The study applied the following criteria; normalised root mean square error (  where O i and P i = observed and predicted value(s), respectively, O i = mean observed data, and x = number of observations. nRMSE defined the developed model's accuracy whilst PBIAS defined the bias provided by the developed model. The error index nRMSE showed the model's performance but did not indicate the degree of over or under-estimation hence the use of the NSE and PBIAS statistical tools in the analysis. The NSE statistic measured the residual variance vs the measured data variance, ranging from −∞ to 1. NSE values between 0.0 and 1.0 are considered acceptable. PBIAS measured the tendency of the simulated data to either under-or over-estimate the observed values. Low magnitudes presented optimal model simulation whilst positive values represented model under-estimation, and negative values represented model over-estimation 33 . The study further employed the prediction efficiency ( P e ). The P e was determined by the R 2 values obtained by regressing the observed and simulated values. A summarised performance rating for the recommended statistics is shown in Table 5.
Scenarios. The developed models were used for situational analysis. The MTI placement depth was maintained at 0.2 m. The operational discharge was varied between 0.27 and 3.19 l h −1 m −1 , which translated to an operating pressure range of 20-160 kPa. The selected discharges used for scenario analysis were selected on the

Results and discussion
Model calibration. The calibration steps for the four wetting geometry models' are outlined below.
Silty clay loam soil. The study followed the outlined steps below to determine the values for the constants A 1 , A 2 , b 1 , and b 2 19 . Recorded volumetric soil water content and observed wetted distance values were used to calibrate the developed soil wetting geometry models. The wetted W and Z were physically measured using grids demarcated onto the transparent plexiglass whilst the volumetric soil water content was measured using the MPS-2 sensors.
Step 1: The dimensionless variables V * , W * , and Z * were estimated using Eqs. (8)- (10) utilising observed values of the requisite variables from the soil bin experiments.
In order to improve the models' accuracy and precision, the calibration step performed iterations on the constants A 1 and A 2 and the resultant equations are shown in Eqs. (22) and (23).   This signified a high-water content in the lateral direction as compared to the vertical direction. This is typical in fine-textured soils (Bouma, 1984). Conversely, the calibration results yielded a scenario where A 1 < A 2 , which subsequently resulted in a W < Z , an observation atypical of finetextured soils. This was because the application times during the experiment promoted the border effect within the confined soil bin.
Sandy soil. The simulation steps for the sandy soil followed similar steps as those described for the silty clay loam soil. The relationships between the dimensionless volumetric soil water content per unit length of MTI ( V * ) and the dimensionless wetted width ( W * ) and wetted depth ( Z * ) are depicted in Fig. 3. The resultant wetted width ( W ) and wetted depth ( Z ) for the soil are shown in Eqs. (21) and (22). W * was plotted against V * , similarly Z * was plotted against V * and the resulting power relationships yielded Eqs. (24) and (25) with R 2 = 0.72 and R 2 = 0.84, respectively (see Fig. 3).
In order to improve the models' accuracy and precision, the calibration step performed iterations on the constants A 1 and A 2 and the resultant equations are shown in Eqs. (28) and (29).
The power indices for the sandy soil, b 1 and b 2 were approximately equal. However, the constant A 1 < A 2 , which resulted in a Z > W , a phenomenon attributed to soil hydraulic characteristics. Gravity forces dominated the soil water movement mechanism in the coarse-textured soil.
It is worth noting that the data in Figs  Considering Eq. (31), a decrease in k by order of magnitude, i.e., migrating to fine-textured soil, yielded a 24% decrease in W and a 119% increase in Z . Likewise, in Eq. (32), a decrease in k by order of magnitude, i.e., migrating to fine-textured soils, resulted in approximately 15% increase in W and an approximately 23% increase in Z.
To assess the sensitivity of the soil wetting geometry with respect to discharge ( q ), the parameters V , D , and k are held constant, and the resultant relationships are outlined by Eqs. (32) and (33).
Doubling q for the silty clay loam soil resulted in a 5% decrease in W and a 27% increase in Z . Similarly, for the sandy soil when q was doubled, there was a 4% increase in W and a 6% increase in Z . Table 6 presents a summarized sensitivity evaluation containing hypotheticals and the resultant wetted horizontal and vertical wetted distances. MTI exhibited an increase in both W and Z in sandy soil whilst it exhibited an increase in Z and a decrease in W for silty clay loam. Regarding the silty clay loam soil, the findings contradict Schwartzman and Zur 15 , who posited that an increase in q results in an increase in W and a decrease in Z of a Gilat loam soil under sub-surface drip irrigation. The difference in behaviour is attributed to the porous nature of MTI, wherein discharge is not of a point source nature.
The sensitivities of W and Z to placement depth D for the respective soils was characterised by Eqs. (34) and (35). An increase in D by a unit magnitude (from 0.2 to 0.3 m) for the silty clay loam soil resulted in an approximately 3% decrease in W and an approximately 53% increase in Z . For the sandy soil, a unit increase in D resulted in a 4% increase in W and a 6% increase in Z . According to Bresler 38 W increases for low k values (fine textured soils) and Z increases by a high magnitude for soils with a high k value (coarse textured soils).
To gauge the sensitivity of W and Z to V, the parameters D , q and k were assumed constant. This yielded Eqs. (36) and (37).

Models' validation.
In order to obtain simulated wetted distances ( W and Z ), an estimation of a range of values for V was made using the sensitivity analysis relationships in Eqs. (36) and (37), and the resultant simulated Z and W for the semi-permeable porous wall line source 2-D flow model were computed. A correlation test based on the R 2 was carried out on a plot of simulated W and Z against observed W and Z (Fig. 4). The correlation coefficients R 2 > 0.75 showed a good agreement between the observed and simulated.
The silty clay loam soil exhibited a wetting pattern on the soil surface, so for that reason, the experimental data from the MPS-2 sensor buried at 0.1 m was excluded. The soil surface wetting phenomenon was observed in both phases of the experiment. A similar observation was made by Kanda, Senzanje 4 for an MTI tubing buried at a depth of 0.2 m. Fan, Huang 39 also made a similar observation and posited that shallow buried depth facilitates upward water movement, a phenomenon observed in fine-textured soils.
The models' evaluation revealed a satisfactory performance. For instance, the silty clay loam models' had a nRMSE of 0.84% and 8.80% for W and Z , respectively, a NSE > 0.5, a PBIAS < ± 25% and an index of agreement ( d ) of 1 and 0.98 for W and Z, respectively (Figs. 5 and 6). The sandy soil exhibited a satisfactory performance as evidenced by a nRMSE of 0.3% and 2.5% for W and Z respectively, a NSE > 0.75, a PBIAS < ± 15% and an index of agreement (d) of 0.6 and 0.3 for W and Z, respectively. The model underestimated the wetted depth ( Z ) for the sandy soil while overestimating the wetted width ( W ). A one-way ANOVA revealed a statistically significant difference ( p < 0.05) in both the observed and simulated wetted W and Z under sandy soil. For the silty clay loam soil, there was no statistically significant difference ( p > 0.05) between observed W and simulated W ; similarly there was no statistically significant difference ( p > 0.05) between observed Z and simulated Z.
Scenario analysis. Silty clay loam. Soil texture influenced water movement ( p < 0.05). Under the silty clay loam soil, the models over-estimated the wetted width ( PBIAS ≤ 18% ), which signified a satisfactory model performance. Other model evaluation metrics (nRMSE ≤ 0.05 , EF < 0.5 ) revealed a good model performance for the wetted (Fig. 7). Interestingly, the model produced a very good agreement between the observed and simulated results at operating pressures of 80-160 kPa. The anticipated result was a good agreement at the manufacturer's design operating pressures of 20-60 kPa. This observation can be potentially attributed to the nature of the Ukulinga soil used for the experiment. The wetted depth (Z) model performed poorly at lower discharges and operating pressures (20-100 kPa). It exhibited a relatively satisfactory performance at operating pressures of (36) Silty clay loam soil    (Fig. 8). Overall, the silty clay loam soil exhibited a pronounced lateral geometry compared to the vertical geometry. The finding concurs with Fan, Huang 39 , who postulated that wetted depth is low in high clay content soils. In addition, a wetting front experiment by Cote, Bristow 11 attributed the fine-textured soil wetting pattern to the dominance of capillarity and matric tension. Saefuddin, Saito 40 observed pronounced radial soil water movement in a silt profile under a multiple outlet ring-shaped emitter. Fan, Huang 39 , in their MTI wetting front experiment, observed a slow lateral water movement after initial wetting of the silt-loam soil after 48 h, much similar to this study timeline.
Sandy soil. Granular soils exhibit a high vertical water movement due to granular pores and dominant gravitational forces; thus, the sandy soil had a greater wetted depth than the silty clay loam soil. The findings concur with a study by Cote, Bristow 11 , who performed a wetting front experiment with a sandy soil under a trickle source emitter. Ghumman, Iqbal 41 attributed the high wetted depths in sandy soils to porosity. Although this study used a line source porous pipe, the findings relate due to the soil hydraulic characteristics factor. Another attribute to high Z in sandy soils is the low macroscopic capillary length in sand which favours gravity as compared to lateral or upward movement. Siyal and Skaggs 42 observed a similar phenomenon in modelling a sandy soil using HYDRUS 2D. The model poorly simulated the wetting distributions under an operating pressure range of 20-160 kPa (Table 7).
Irrigation implication. The developed soil wetting geometry models can be adopted to estimate soil wetting geometries for particular soils in question, thus ensuring optimal placement depth and spacing of MTI laterals. Knowledge of wetting geometries can potentially aid irrigators in adopting optimal lateral placement depth and spacing. The wetting depth models can inform irrigators on irrigation application times, thus availing optimal volume to the root zone and minimizing deep percolation and soil water loss due to evaporation. For fine-textured soil, shallow buried depth avails water to the soil surface that will be lost due to soil evaporation. The lateral water front for fine-textured soils expanded faster than that of coarse-textured soils; hence lateral  www.nature.com/scientificreports/ placement under fine texture soils should be strategically placed to promote an optimal overlap between row crops. For coarse texture, close lateral placement will be required to create an optimal wetting front overlap. (21) and (24)-(25) may require testing in different geographical locations and assess their universal suitability. The study was conducted for 20 h under a silty clay loam and 96 h for the sandy soil. The soil bin's lateral dimension (width) limited the testing times for silty clay loam. The implication was the influence of border effects on soil water movement if the experimental times went beyond 20 h. Likewise, for the macroscopic sandy soil, the experimental times were limited to 96 h because of the soil bin's depth restrictions. Models' development was also limited to the following constant inputs; placement depth ( D = 0.2 m) and discharge ( q = 2.39 l h −1 m −1 ). The wetted depth Z is not entirely an independent variable as factors such as crop-specific root water uptake influence the vertical soil-water movement.

Conclusions and recommendations
The study adopted the Buckingham π theorem or dimensional analysis to develop, calibrate and validate models' that simulate soil wetting geometries for MTI as a function of soil hydraulic conductivity ( k ), placement depth ( D ), emitter discharge ( q ), and soil water content ( V ). Soil texture significantly affected the wetting geometry under MTI. The silty clay loam models satisfactorily simulated the wetted width or lateral soil moisture movement. It, however, failed to simulate the wetted depth at an operating pressure range of 20-100 kPa, which translated to an operating discharge of 0.27-1.18 l h −1 m −1 . It can be concluded that the models are prescriptive, i.e., they should be used for the soils and conditions that they were developed under.
The study also noted, judging from the wetting pattern of the fine-textured soil, there is potential in MTI to provide plants with water with minimized deep percolation losses. The study was done in a soil bin on bare homogenous soil. The researchers recommend the study be carried under field conditions for both cropped and un-cropped soils and test the suitability of the developed models. Furthermore, the experiment was carried on dry soil. Thus an investigation on the soil wetted pattern under an initially moist soil should be carried out and compared to the current study ( Supplementary Information 1).