A Computational Method for Modeling Spatiotemporal Variability of Hydrodynamic Properties in Sandy Soil Under Drainage and Recharge

This article proposes an analytical strategy that combines X-ray Computed Tomography (CT) and Continuous Wavelet Transform (CWT) analysis as an alternative solution to long-term experiments that seek to investigate spatiotemporal variations in soil hydraulic properties induced by drainage and recharge cycles. We conducted CT scanning on 100-cm-high column filled with two types of sandy soil in a laboratory environment to simulate, over the period of a month, the equivalent of nearly 40 years of drainage/recharge cycles akin to agricultural fields adopting subirrigation as water management practices. We also monitored soil matric potential, water inflow and outflow, as well as movement of tracers. This later consists in zirconium oxide (ZrO2) that we added to the top 20 cm of each soil column. The results revealed that drainage and recharge cycles greatly affect the evolution of soil hydraulic properties at different locations along the soil profile by reducing drainage and capillary capacities. The approach also allowed us to identify each periodic component of drainage and recharge cycles, and thereby calculate the periodic drift over time. The proposed method can be applied to predict soil evolution according to soil texture, drainage system design and water management, thereby offering a potential basis for proposing mitigation measures related to soil hydrodynamics. It may find its application in agricultural farms adopting subirrigation and surface (e.g., drip) irrigation approaches and, in mining and civil engineering.


INTRODUCTION
Subsurface drainage systems, a commonly used agricultural water management practice, can considerably influence the formation and evolution of soils (1,2). Subsurface drainage systems alter the frequency of water table fluctuations (drainage and recharge cycles) and soil physicochemical properties (3,4) including soil pores size distribution (5,6). The repeated drainage and recharge cycles reduce percolation, induce soil shrinkage, and decrease soil bulk density, as observed in rice paddies (7), soybean fields (2), and for a range of soil types and land use conditions (8). These effects occur over both space and time (9,10) and intensify alterations of soils under anthropogenic activities (11). An example of such anthropogenic soils is observed in cranberry fields which consist of natural or constructed sandy and peat soils.
The spatiotemporal variability in soil profiles induced by drainage systems is capable of influencing yield and irrigation water use. Gumiere et al. (12) found that: (a) areas of low crop yield are spatially correlated with soil horizons of low saturated hydraulic conductivity, and (b) accounting for spatial variability of soil hydraulic properties holds the potential to reduce irrigation water use by up to 75%, in a highly drained sandy soil under cranberry cultivation. Therefore, the prediction of spatiotemporal variations in soil hydraulic properties under drainage conditions, especially in cranberry fields with subirrigation, is of prime importance.
Assessing spatiotemporal variation in soil hydraulic properties requires laboratory analysis of soil samples in addition to in situ observations (13). However, core sampling may alter soil properties, such as bulk density and pore size distribution, leading to differences between in situ and laboratory measurements of soil water retention and unsaturated hydraulic conductivity (14). Similarly, soil properties may vary rapidly following tillage, a process that standard soil sampling techniques may fail to capture (15). Furthermore, in situ experiments often require extensive and continuous field instrumentation as well as periodical instruments removal to allow agricultural practices (16,17). The time-consuming and labor-intensive nature of these field operations suggest the need for more convenient and, more importantly, non-destructive techniques to examine soil structures.
Various non-destructive imaging techniques have been proposed as solutions to the abovementionned problems with in situ soil structure analysis. The devices that implement these techniques work by producing signals that contain information on the targeted soil properties that are subsequently inferred using a statistical model. Electrical resistivity tomography can characterize active zones of water uptake (18). Spectral electrical impedance tomography and electrical impedance spectroscopy have been used to characterize and monitor crop root systems (19). Ground-penetrating radar (GPR) can measure soil water content (20). A thorough review of soil imaging methods and their corresponding measurable soil properties can be found in (21).
Among non-destructive techniques, X-ray computed tomography (CT) is acknowledged as a technique to observe and measure soil structure (i.e., shape or type, size) and soil pores characteristics (i.e., quantity, size, location, and shape), and soil rupture resistance (i.e., plasticity, toughness, stickiness, penetration resistance and excavation difficulties). CT has been used for high-resolution monitoring of time-dependent changes in soil properties, such as bulk density, porosity, volumetric water content, and solute transport parameters (22), aggregate characteristics (23), and water characteristic curves (24). A detailed account of X-ray CT technology, and its application in soil science can be found in (25,26), respectively.
We hypothesized that a combination of imaging (e.g., Xray CT scanner) and non-imaging in situ sensors (e.g., soil matric potential) can yield time series of the data required to quantify changes in soil drainage and capillary capacities. However, the non-stationary nature of these time series calls for statistical models such as Fourier transforms (FT) and wavelet transforms (WT) (27). Prior studies have suggested that WT are more effective than FT in analyzing non-stationary data as they provide a time-scale representation of the data (27). Therefore, as with the spectral analysis performed by (5), WT can decompose a time series, such as soil matric potential, into its contributing components, i.e., periods of drainage and recharge. Hence, it enables a complete representation of the timings between drainage and recharge events with different frequencies of occurrence.
The main objective of this study was to define a method for analyzing long-term spatiotemporal variability in hydraulic properties of a sandy soil under repeated drainage and recharge cycles without waiting several years before obtaining the results. To this end, soil column experiments were performed to: (a) simulate nearly 40 years of water-table depth variations in sandy soil samples obtained from a cranberry field in Quebec (Canada), through down-scaled and accelerated, drainage and recharge cycles, and (b) assess the soil hydrodynamic response. The proposed strategy aims at anticipating changes in soil hydraulic properties that may occur following the construction of cranberry fields equipped with subsurface drainage systems. It represents an alternative to long-term field experiments and traditional methods for assessing potential spatiotemporal variations in soil hydraulic properties.

Soil Columns Preparation
At the INRS CT Scanning for Civil Engineering and Natural Resources Laboratory 1 (Quebec, Canada), we set up an experiment consisting in four transparent plexiglass soil columns (height: 100 cm, diameter: 15.24 cm- Figure 1A) packed with sand. The percentage of silt and clay were negligible. The columns contained either a fine or medium-coarse sandy soil or the two grades in 20-80% proportions. The soil was added into the columns in 2-cm-thick layers at a uniform bulk density of 1.6 g/cm 3 . The upper 20 cm of each column was a mixed with zirconium oxide (ZrO 2 ) tracer (5 g of ZrO 2 per 100 g of sand). The ZrO 2 tracer was used to simulate the movement of finer, but heavier particles compared to sand particles, and consisted of spherical colloidal particles between 0.01 and 10 µm in size, with a median particle size (d50) of 1.5 µm and was considered inert. The ZrO 2 particle size distribution is shown in Figure 1B.

Experimental Procedure
We performed experiments involving four soil profiles: • Type-1 composed of fine sand only (d50 = 174 µm) and named FS hereafter,  Each experiment started with a recharge event that raised the water level in the column to a depth of 35 cm below the soil surface, followed by a drainage event that lowered the water level to a depth of 55 cm. To monitor the water table fluctuations, we installed 10 tensiometers equipped with pressure transducer (PX-26, Omega Engineering Inc., CT, USA) at 10 cm intervals between depths of 5 and 95 cm from the top of each column ( Figure 1A). The tensiometers were connected to the control unit; recording average soil matric potential every 5 min.
Recharge was automatically activated as soon as the tensiometer installed at 55 cm recorded a value of soil matric potential below zero, indicating that the water table had reached that depth. Similarly, recharge was stopped, and drainage activated as soon as the tensiometer installed at 35 cm recorded a value below zero. The number of recharge and drainage cycles, as well as the start and stop time of drainage and recharge events, were automatically recorded by a CR10 control unit (Campbell Scientific, Logan, UT, USA). The simulation duration, which is 38.92 years, was determined by measuring soil matric potential time series during a growing season and deducing that a year corresponds to approximately 75 drainage/recharge cycles. We therefore calculated the number of years corresponding to each simulation by dividing the cumulative number of cycles by 75. This approximation is conservative because we considered only the number of cycles during the growing season, which corresponds to the periods of artificial drainage/recharge events in cranberry fields.

Soil Matric Potential Time Series Analysis Using Continuous Wavelet Transform
The Continuous Wavelet Transform (CWT) of a time series vector y is defined as the convolution of y m with a scaled and translated version of the wavelet function ψ: where the amplitudes of wavelet representations (given by the complex conjugate ψ * ) of the original time series y are scaled and translated along the continuous wavelet scale (s) and discrete localized time index (m). The function ψ must be localized in both time and frequency and have a zero mean to qualify as a wavelet function (28). Here, we used the Morlet wavelet (29), a complex non-orthogonal plane wave (i.e., a wave with equally spaced maxima) modulated by a Gaussian function: where w 0 is the non-dimensional frequency and η is a nondimensional time parameter. We used w 0 = 6 to satisfy the zeromean condition. CWT analysis is facilitated by a wavelet power spectrum which can be obtained by normalizing the absolute expected value of the Morlet wavelet transform (W m (s)) with respect to the series variance σ 2 over a range of scales: In using the normalized power spectrum to perform CWT analyses, the longest period for which a trend can be detected depends on the length of the time series. In this study, we applied the CWT to soil matric potential time series recorded by the tensiometer installed at 15 cm below the soil surface. This depth was selected because it is just above the interface between the two soil horizons used in the FS/MCS and MCS/FS columns and corresponds to the bottom of the cranberry root zone (30). Frontiers in Soil Science | www.frontiersin.org Our objective in using the CWT algorithm was to quantify the duration between occurrences of two successive drainage and recharge events. This analysis was performed using the library WaveletComp (31) of the R programming language (32).

CT Scanning
The experiments were performed using a Somatum Volume Access CT scanner (Siemens, Munich, Germany). Two energy levels were used, 140 and 120 keV, and the spatial resolution (voxel) was 0.06 × 0.01 × 0.01 cm 3 (x, y, z). The CT images were recorded in DICOM V5 format and analyzed using the R programming language (32) with the oro.dicom library (33). The spatiotemporal variations in soil properties were monitored on a weekly basis over a period of 7 weeks. Before each scan, the soil columns were completely saturated through irrigation activation.

Voxel Porosity
Changes in the soil porosity and tracer (ZrO 2 ) concentration can be obtained from the X-ray CT scan data using the following equations (34): C s (z) = HU 1 (z)HU Zr2 − HU 2 (z)HU Zr1 HU Zr1 HU m2 − HU Zr2 HU m1 (5) where C Z r is the tracer concentration at depth (z); C s is the solid phase concentration (assumed to be equal to the concentration of quartz) at depth (z); φ is the porosity as a function of depth; HU m1 and HU m2 are the absorption coefficients of quartz at energy levels of 140 and 120 keV, respectively, which correspond to 1,798 and 1,886 HU (Hounsfield Units); HU Zr1 and HU Zr2 are the absorption coefficients of ZrO 2 at energy levels of 140 and 120 keV, respectively, which correspond to 14,170 and 16,650 HU; and HU 1 and HU 2 are the mean absorption coefficients (image slices from the X-ray CT scans of the soil columns) of saturated soil at energy levels of 140 and 120 keV, respectively.

Unsaturated Hydraulic Conductivity of the Soil
We used the Carnahan-Starling approximation of the near-surface complementary cumulative density function (cdf) for voids [e vi (δ, z); (35,36)] to derive the water Frontiers in Soil Science | www.frontiersin.org retention and unsaturated hydraulic conductivity curves of the soil: where δ is the particle radius (µm); η (z) is the dimensionless density, which is calculated as 1φ(z); S i is the surface area ratio; m 1i is the first moment of the dual lognormal particle size distribution used in this study, and a 0i , a 1i , and a 2i are coefficients corresponding to the moments of the distribution.

Measured and Reconstructed Time Series of Soil Matric Potential
The soil matric potential measured at 15 cm with respect to the soil surface for all soil columns used in this study are presented in the upper panels of The results also demonstrate the differences in drainage behavior between homogeneous and heterogeneous soil profiles as well as between coarse-over-fine and fine-over-coarse soil layers as reported by (37,38). The initial oscillation, followed by the stabilization in the time series as shown in Figures 2, 5, reflects rapid changes in soil consolidation occurring in cranberry fields immediately after construction. We were able to reconstruct the original time series from the wavelet transform with a coefficient of determination (R 2 ) close to 0.9; indicating that the wavelet transform was able to adequately reproduce the time evolution of the main oscillating components of the time series.

Wavelet Power Spectra
The resulting wavelet power spectra of the experimental drainage and recharge cycles are presented in the lower panels of  Frontiers in Soil Science | www.frontiersin.org periods, i.e., drainage. They reflect the duration of drainage and recharge events.
For all soil profiles, the spectrograms generally show an increase in the periods of the red and blue bands over time, indicating an increase in duration between successive drainage and recharge events. This may be due to the reduction in porosity and to the accumulation of small particles in the pores occurring during water table fluctuations within the soil profiles and they are the governing processes responsible for the increase in duration of drainage and recharge cycles within the soil profiles. Due to the spatial nature of drainage and recharge processes, this increase in duration suggests a spatiotemporal shift toward lower values of soil hydraulic properties; an indication that the proposed method bears the potential to simulate spatiotemporal variations in soil hydraulic properties.
All columns, except the MCS/FS column, present high period shifting. The periodic drift for FS/MCS profile did not change much although FS/MCS soil profiles are prone to the suffusion phenomenon which is an internal erosion of the fine particles through the voids of the coarse ones due to seepage flow.
For the FS and MCS soil profiles, a halt in the cycle duration increase is observed at different points in time (lower panels of Figures 2, 5). This is because it became impossible for the system to reach the minimum and maximum matric potential thresholds beyond these points in time. This implies considerable changes in soil hydraulic properties presumably due to the reduced capillary capacity caused by the particle size discontinuity in these columns, preventing the water table from reaching a matric potential triggering drainage.
Recharge and drainage periodicity, as a function of simulation time, is presented in Figure 6. We observed an increase in the periods between successive drainage and recharge events except for the Type-2 profile (Figure 6B). This behavior may be due to the dynamic characteristics of the soil during its evolution, i.e., the migration of small particles and the reorganization of coarse particles in the soil matrix during recharge and drainage.

Spatiotemporal Variability in Total Porosity
We present the spatiotemporal variability (weekly) in total porosity (φ) and ZrO 2 concentration for the soil columns in Figure 7. We observed a reduction in φ at the bottom of the sand-ZrO 2 layer mix, i.e., at a depth of approximately 20 cm, for MCS/FS and FS/MCS profiles only, but no noticeable reduction below this layer for the other soil profiles. This was probably due to the less intense soil packing at the interface between the two heterogeneous (fine-coarse) layers that composed these two profiles. However, there is a slightly different pattern in  total porosity change below the sand-ZrO 2 layer between FS and MCS profiles. For instance, there is an increase in φ below the interface for the FS profile. These results can be explained by the downward migration of small particles, especially ZrO 2 .
The movement of small particles may have destabilized the soil structure at approximately 20 cm and increased φ due to the reorganization of coarse particles. This increase in φ is consistent with the increases observed in the size and number of coarse pores during wetting-drying cycles reported by (5).

Colloidal Particle Transport
We observed the migration of small particles in all four soil columns. We noticed that the ZrO 2 concentration increased in the upper part of the water table fluctuation zone. It is important to note that before each scan, the soil columns were completely saturated by a recharge event, complemented when necessary with irrigation. Therefore, the upward water flow may have caused upward movement of small ZrO 2 particles. Indeed, the ZrO 2 particle size distribution was not uniform (Figure 1) and the fine size fraction used in the experiments can be very mobile.
Below the surface horizon (0-20 cm), the ZrO 2 concentration profiles were indicative of hyperexponential deposition. This phenomenon is caused by straining and mechanical filtration, which primarily occurs near the outlet source of colloids and at the smallest pores network (39). Considerable downward migration of ZrO 2 was observed, and this migration caused the observed increase in the corresponding ZrO 2 concentration at depths from approximately 20 cm to the bottom of the column for MCS/FS and MCS profiles. The increase in concentration indicates accumulation that may be attributed to clogging, straining, and attachment to solid soil faces and to the air/water/solid interface (40)(41)(42).
Our results corroborate with those of (43,44), who observed greater compaction under drainage conditions in coarse sand compared to fine sand, due mostly to a low capillary force caused by the presence of large pores. They are also consistent with those of (45), who observed a reduction in the bulk density measured by gamma-ray CT scans during experiments involving wettingdrying cycles. Ma et al. (46) also observed that wetting-drying cycles affect the porosity, mainly by reducing macropores with diameters greater than 100 µm.

Changes in Soil Hydraulic Properties
For each of the four soil columns, Figure 8 shows the unsaturated hydraulic conductivity curves for drainage (top set of curves) and recharge (bottom set of curves) at the locations that experienced the greatest porosity reductions in the soil profiles. At low soil matric potentials and near-saturation conditions, unsaturated hydraulic conductivity decreased rapidly with increasing matric potential, which is typical of sand. These large changes in unsaturated hydraulic conductivity of the soil columns contributed to modifying the drainage and capillary capacities of the soil, as evidenced by the drops in the periods of drainage and recharge cycles given by the CWT analysis (see for example Figure 3). Furthermore, the pair of curve sets in Figure 8 depicts the effect of the hysteresis phenomenon on the hydraulic behavior of the soil profiles. The results show that hysteresis does not affect the porosity for the studied soil profiles, but is a phenomenon solely dependent on soil moisture and hydraulic conductivity.
In real-world situations, different factors, including management practices, biological factors, rainfall, soil consolidation or sedimentation, soil disaggregation, drying and wetting processes, initial soil moisture, and erosional and depositional processes, have been identified to explain changes in soil hydraulic conductivity (2). The soil is relatively loose following tillage and irrigation, thus, the first infiltration event for various textured soils is another factor that can also induce strong variations in hydraulic properties based on the soil type and irrigation method (5). The soil gradually compacts during subsequent wetting and drying cycles. For instance, Zeng et al. (47) observed relative rates of change of 0.57 and 0.84 in the saturated hydraulic conductivity (Ksat) of sandy loam under furrow irrigation (2,250-2,700 mm of water) and drip irrigation (675-787 mm of water) after 15 irrigation events, respectively. Our study showed similar results, i.e., a reduction of Ksat over time for all columns except for the Type-2 profile. The water pressure boundary of 9 cm during irrigation and complete saturation is similar to the conditions found in paddy soil (48). In particular, Yoshida et al. (48) and Zhang et al. (49) studied paddy soils and found that the soil drainage capacity is reduced under submerged conditions associated with flooding as a result of diminished hydraulic conductivity and that this reduced drainage is caused mainly by dispersion, flocculation, migration and reorganization of soil particles.
The proposed method did not consider the effects of varying temperatures and pressures and the influence of evaporation; therefore, further work should extend the method to consider these conditions. Also, it should be extended to other soil types and management practices, in addition to research in other types of irrigation methods such as drip irrigation.

Practicalities of the Proposed Method
The results presented in this article further demonstrate the relevance of CT and CWT in anticipating the spatial variability of anthropogenic-soil properties. However, the typical cost of running the experiments may be a prohibiting factor to widespread application and replication of the proposed method. Fortunately, nowadays, X-ray CT has become more commonplace in the earth sciences for imaging soil samples (25,26,50). For instance, recent developments in X-ray microtomography have led to many commercially available and miniaturized desktop CT systems, making CT and even micro-CT (higher spatial resolution typically around 1 µm) accessible to a large number of researchers (51,52). Furthermore, industrial X-ray CT scanners are typically custom-built (25), which makes it possible to, for example, repurpose a scanning electron microscope into an X-ray microscope-microtomography to generate X-rays for tomographic reconstruction in the micrometer range (51). Such a solution, although somewhat limited by the recording time, may be a very affordable alternative to adequately replicate our experiments. Also, inexpensive and compact CT units may not be able to accommodate the full length of the plexiglass columns used in this study because of their smaller field of view. This limitation may be overcome by applying an appropriate tomographic images stitching strategy as reviewed by (53). A recent example demonstrating the application of such strategies, termed Tomosaic, is presented by (54).
This study has been conducted in a laboratory environment and has not considered a number of real-world factors. These factors include microbial and root effects as well as field interventions that add sand layer to foster rooting, which are hard to simulate. This limitation suggests the need for fieldbased experimental verification of the computational method over a longer period of time. This is all the more important to verify the effects of raising the water table to much higher levels and the associated speed would have on colloidal particle movement, although the conclusion of this study is consistent with the amount of compaction observed in the field.

CONCLUSIONS
We proposed an analysis strategy based on X-ray CT for predicting soil water retention and hydraulic conductivity, particle size distribution, and the CWT of the soil matric potential. The main advantages of the proposed method are its ability to: (i) calculate decreases in porosity and unsaturated hydraulic conductivity for drainage and recharge cycles and small particle concentrations at a high resolution (600 µm); (ii) rapidly conduct analyses relative to traditional long-term field experiments; (iii) simulate years of drainage and recharge cycles; and (iv) identify each component of individual drainage and recharge cycles and calculate the periodic drift in time. The proposed method presents an alternative solution to long-term experiments in the field and traditional methods for determining spatiotemporal variations in soil hydraulic properties. It may allow the evolutionary characteristics of soils to be predicted according to the vertical distribution of particle sizes, drainage system design and water management practices implemented, thereby allowing variations in soil hydrodynamic parameters to be anticipated and possibly controlled.
The method is equally applicable to non-agricultural fields subjected to naturally occurring drainage and recharge cycles. Furthermore, we anticipate a potential application in dripirrigated fields as the high frequency of the water application under drip irrigation influences soil properties as well. The method is also relevant to the mining industry as wetting and drying cycles influence, for example, the mechanical strength characteristics of coal samples and the microstructure of rock mass. The civil engineering industry is another potential field of application as the strength of waste materials such as clay and fly ash used for developing sustainable construction materials may be impacted by wetting and drying cycles.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
YP, JL, and SG were responsible for the planning and design of the research work. YP and SG carried out the experiments and data analysis. PC participated in the drafting of the manuscript and interpretation of the results. JL, JG, AR, JC, and TG reviewed and proofread the article. All authors contributed to the article and approved the submitted version.