Continuous monitoring of a soil aquifer treatment system's physico-chemical conditions to optimize operational performance

Soil aquifer treatment (SAT) is a tertiary process for wastewater treatment, where the wastewater infiltrates through a thick vadose zone for purification and storage in the underneath aquifer. SAT infiltration basins are typically flooded intermittently, while maintaining a fixed ratio between the wetting and the drying stages. However, infiltration basins exhibit different physical and chemical properties, limiting the generalization of SAT operation to attain optimal efficiency. Since frequent sampling of the soil pore water to verify the SAT's biodegradation efficiency can be arduous, continuous monitoring of the SAT vadose zone's physico-chemical conditions is required. In this study, redox potential (E_h) was continuously monitored, together with other variables, such as volumetric water content (θ), soil temperature, and gaseous oxygen (O₂), at multiple depths of a SAT vadose zone throughout the year and while the system was constrained to different operational modes. Hydrological models were calibrated and validated to water content observations, and they illustrated the seasonal changes in water infiltration. Furthermore, it was shown that, under long wetting stages during winter, there was a reduction in the SAT's drainage capabilities. The E_h observations, under long wetting stages, demonstrated larger variability and very negative values as ambient temperature increased. Assembling the daily E_h observations illustrated that a wetting stage should cease after about 30 h, once suboxic conditions are established. A drying stage's optimal duration should be 36 h, according to the E_h and O₂ observations during summer and winter. Ultimately, the study shows that the length of wetting and drying stages should be defined separately, rather than by adhering to the wetting $/$ drying ratio.

1 Introduction

Worldwide water scarcity has motivated the development of alternative water resources such as the reuse of treated wastewater. Soil aquifer treatment (SAT) is commonly implemented to further improve the recovered water's quality and remove the majority of suspended matter, microorganisms, viruses, and organic and inorganic constituents (Dillon, 2005; Goren et al., 2014; Massmann et al., 2006; Schmidt et al., 2011; Tsangaratos et al., 2017). In SAT systems, the treated wastewater is recharged to the underlying aquifer by surface spreading over infiltration basins. The wastewater is purified mainly through the physical and biochemical processes that occur during water passage through the vadose zone (Dillon, 2005; Elkayam et al., 2015). Although SAT systems have been used for decades (Grinshpan et al., 2021; Bouwer, 2002), the ability to estimate and predict a SAT system's performance is still challenging, and the optimal SAT operation is still under investigation (Ben Moshe et al., 2020; Sharma and Kennedy, 2017).

A major uncertainty in SAT systems concerns the vadose zone processes that play a central role in determining the quality of the water that recharges the aquifer (Elkayam et al., 2015). The chemistry of the percolating wastewater changes due to a combination of several biogeochemical processes, such as organic matter biodegradation, nitrification, sorption, cation exchange, etc. (Amy and Drewes, 2007; Díaz-Cruz and Barceló, 2008; Goren et al., 2014; Miller et al., 2006; Tufenkji et al., 2003). Most of the organic matter is removed by biodegradation (i.e. microbial activity) within the upper two metres of the vadose zone (Drewes, 2009). Nevertheless, the microbial activity is greatly affected by the soil water content, which frequently changes in SAT. Generally, a major challenge in SAT systems is to facilitate the intrusion of O₂, primarily in the gaseous phase, and to enrich the active subsurface with O₂ (Ben Moshe et al., 2020; Massmann et al., 2006).

A consequence of the perturbation in the O₂ supply to SAT is expressed in changes in the redox conditions (Mächler et al., 2013; Rezanezhad et al., 2014). Redox potential or oxidation reduction potential is a quantitative measure of electron availability, i.e. the tendency of the system to receive or donate electrons (Hinchey and Schaffner, 2005). Substantial changes in SAT systems' redox conditions might lead to the release of undesirable metals, such as Fe²⁺ and Mn²⁺ (Goren et al., 2012), and affect the degradation rates of pesticides and pharmaceutical substances (Massmann et al., 2006). Additionally, previous studies have illustrated the possible degradation of groundwater quality due to the emergence of contaminants that leach from the SAT vadose zone under reducing conditions (Asano and Cotruvo, 2004; Massmann et al., 2006; Oren et al., 2007; Sharma and Kennedy, 2017). Redox processes are associated with the degradation of organic matter by terminal electron acceptors or redox couples, such as O₂ $/$ H₂O, NO₃ $/$ N₂, MnO₂ $/$ Mn²⁺, Fe³⁺ $/$ Fe²⁺, and SO₄ $/$ H₂S, in sequential order from the highest energy yield downwards (Berner, 1981; Froelich et al., 1979; Christensen et al., 2000). The transition between redox conditions is determined by the presence and availability of these electron acceptors/donors. Once the strongest oxidizing species (O₂) is depleted, the next strongest oxidizing species is used (NO₃), and so on. The alternation between oxic (> 400 mV), suboxic (between 400 and 100 mV), and anoxic (100 mV >) conditions in the vadose zone depends on the availability of the oxidized species (Reddy et al., 1998). In addition, studies have reported on the seasonal (temperature changes) effects on redox conditions, which were attributed to the increase in dissolved oxygen concentrations at low temperatures (Massmann et al., 2006) and the greater microbial activity (i.e. higher O₂ consumption) at higher temperatures (Greskowiak et al., 2006; Kirschbaum, 1995).

An important operational aspect of a SAT system is the intermittent application of the effluents (Sattar, 2016; Sallwey et al., 2020). After the infiltration basin is flooded with wastewater, a drying period is implemented to sustain the SAT's infiltration capacity and biochemical capabilities (Sharma and Kennedy, 2017). The wetting and drying stages, which can be expressed by the wet $/$ dry ratio parameter, have a critical impact on the removal rates of dissolved organic carbon, total nitrogen, and pathogens (Ben Moshe et al., 2020; Morrison et al., 2020; Sharma and Kennedy, 2017). Although the wet $/$ dry ratio can vary depending on location and wastewater quality, it is well accepted that it should be below 1.0 (Sattar, 2016; Sharma and Kennedy, 2017). Nevertheless, infiltration basins behave differently with regard to infiltration rates and clogging. Thus, in many cases, the SAT operational efficiency is limited to the personal experience of the operators and their knowledge of the specific infiltration basin (Sharma and Kennedy, 2017). Note, however, that several studies (e.g. Ben Moshe et al., 2021) suggest that it is not the wet $/$ dry ratio that should be considered but specific times for wetting and drying.

The oxidation reduction potential (E_h), together with chemical and physical parameters, such as water content, soil temperature, O₂ concentration, etc., can be continuously monitored by installing the relevant sensors. Previous studies have implemented the E_h sensor and successfully described, with high temporal resolution, the subsurface chemical conditions in various environments, such as wetlands, the groundwater (or the capillary fringe), aquifers, etc. (Wallace et al., 2019; McMahon and Chapelle, 2008; Shenker et al., 2005; Silver et al., 2018; Rezanezhad et al., 2014). To improve SAT system performances, the link between the wetting and drying stages and the subsequent redox conditions developed in the subsurface should be established. Thus, in situ monitoring can improve SAT management performance and reduce the subjectivity of the operator. The objective of this study was to examine the temporal variability in redox potential and the way it is affected by changes in volumetric water content, gaseous O₂, and climate imposed by different operational modes of wetting and drying stages. Furthermore, calibrated and validated hydrological models were used to explore the behaviour of water fluxes under different operational modes and seasonal temperature changes. Finally, the optimal lengths of a drying stage and a wetting stage were determined, following the in situ observations.

2 Methods

2.1 Study sites

The Dan Region Reclamation Project (Shafdan) reclaims about 125 × 10⁶ m³ of effluent annually from the Tel Aviv metropolitan area in Israel. The treatment of effluents occurs in two stages. The first stage involves mechanical biological treatment, which is based on activated sludge, while, in the second stage, the treated water (a secondary effluent) is delivered to infiltration basins, as part of the SAT system, to further improve water quality. The six infiltration basin sites, covering a total area of 1053 km², are located in central Israel, overlying the coastal aquifer (Fig. 1). Each basin is divided into several spreading ponds, about 1500 m² each, which are alternately flooded. The vadose zone that underlies the basins is mostly composed of sand, sandy loam soil, and calcareous sandstone layers. Typically, the ponds are flooded for 1 to 2 d (max hydraulic head of about 50 cm), followed by 2 to 6 d of drainage and soil surface drying. The wetting and drying stages are controlled by the ponds' flooding order, the availability of effluent, and the drying period, which is suggested to be at least 24 h (Icekson-Tal et al., 2003). The basin surface is ploughed on a regular basis to break up the developed biocrust and to prevent clogging (see Negev et al., 2020, for details).

2.2 Study operation

In this research, two in situ measurement stations were installed in an infiltration pond during 2018 (pond 4103 in the Yavne 1 cluster; Fig. 1). Each station was equipped with several sets of sensors at 25, 50, 75, 100, and 150 cm depth, including time domain transmittance (TDT) probes (Acclima Inc., Idaho, USA), copper–constantan thermocouples (Omega Engineering Inc., CT, USA), oxidation reduction potential (ORP) electrodes (ELH016; Van London Co., Houston, TX, USA), and O₂ percentage probes (SO-421 sensor; Apogee Instruments, Utah, USA). Data were collected at prescribed intervals and logged on a CR1000 data logger (Campbell Scientific, Inc., Logan, UT, USA). In addition, suction cups were installed at similar depths. In station 1 (Fig. 1), the data consisted of volumetric water content (θ), soil temperature (T), and ORP (E_h) time series, which were continuously measured every 20 min between 28 July 2020 and 10 February 2021 (a total of 14 185 values, over 197 d, for each variable). The data were obtained at 25, 50, and 100 cm depths. In station 2 (Fig. 1), θ, T, gaseous oxygen (O₂), and E_h were measured every 20 min between 8 May 2019 and 20 July 2020. There were about 60 d in which data were not collected in station 2 due to technical issues. The data from station 2 contained 27 222 points of θ, 29 394 points of O₂, 30 414 points of E_h, and 26 730 points of T measurements. In station 2, the data were collected at 25, 50, 75, and 100 cm depths.

The water quality characteristics of the secondary effluent that flooded the Yavne 1 basin are presented in the Supplement (Fig. S1). Note that the quality parameter concentrations conform to the updated Inbar regulations (Inbar, 2007) and the findings of a previous study that surveyed numerous wastewater storage and treatment reservoirs across the country (Kfir et al., 2012). To determine the soil physical properties, undisturbed soil cores were sampled at different depths. Subsequently, flow experiments were carried out to calculate the saturated hydraulic conductivity (Ks), based on Darcy's law (Fig. S2). Additionally, particle size distribution (PSD) analyses are presented in Fig. S3. The PSD results indicate that the SAT vadose zone is homogeneous.

Figure 1 The location of the investigated site – the Yavne 1 infiltration basin of the Shafdan. In the magnified view of the investigated pond, the yellow circles represent the locations of the measurement stations (© Google Earth).

2.3 Hydrological model and gaseous oxygen dynamics

The calculations for water and oxygen fluxes in the SAT vadose zone are calculated differently for the ponding stage and for the stage, where there is no ponding (water) on the soil surface. For the ponding stage, Ganot et al. (2017) showed that the infiltration rates in managed aquifer recharge systems can be predicted reasonably well by simple analytical models. In this study, the Green–Ampt equation for infiltration into a flooded soil was implemented to calculate the water flux as follows (Bouwer, 2002):

$$\begin{array}{}\text{(1)}& q=\mathrm{Ks}\times {\displaystyle \frac{\left(L+d-{\mathit{\psi }}^{*}\right)}{L}},\end{array}$$

where q (L T⁻¹) is the water flux, d (L) is the ponding depth, L (L) is the thickness of the saturated vadose zone, and ψ^* (L) is the negative pressure head at the wetting front. Note that L is assumed here to be constant, the subsurface is assumed to be homogeneous, and θ is assumed to vary with time only. As the wetting front progresses, the gradient approaches a value of unity, and the infiltration rate becomes equal to the hydraulic conductivity of the wetted zone. Once the water ponding ceases, the water drainage is set equal to the unsaturated hydraulic conductivity, described with an exponential form as follows (Guswa et al., 2002):

$$\begin{array}{}\text{(2)}& D\left(\mathit{\theta }\right)=\mathrm{Ks}{\displaystyle \frac{e^{\mathit{\beta }\left(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{fc}}\right)}-\mathrm{1}}{e^{\mathit{\beta }\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{fc}}\right)}-\mathrm{1}}},\end{array}$$

where Ks (L T⁻¹) is the saturated hydraulic conductivity, β is a parameter of the soil, θ (L³ L⁻³) is the volumetric water content, θ_fc (L³ L⁻³) is the water content at field capacity, and θ_s (L³ L⁻³) is the saturated water content. Furthermore, the effect of temperature changes on the soil hydraulic conductivity is implemented through the change in viscosity as follows (Lin et al., 2003):

$$\begin{array}{}\text{(3)}& {\mathrm{Ks}}_T={\mathrm{Ks}}_{\mathrm{25}}{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{25}}}{{\mathit{\mu }}_T}},\end{array}$$

where Ks_T and Ks₂₅ are soil hydraulic conductivity values at temperature T and 25 ^∘C, respectively, and μ_T and μ₂₅ are the dynamic viscosity of water (M L⁻¹ T⁻¹) values at T and 25 ^∘C, respectively.

An inverse problem was set to find an optimal combination of Ks and β parameters that minimizes the following objective function:

$$\begin{array}{}\text{(4)}& \mathrm{\Phi }\left(b\right)=\sum _{i=\mathrm{1}}^N{\left[\mathit{\theta }\left(t_i\right)-\mathit{\theta }\left(t_i,b\right)\right]}^{\mathrm{2}},\end{array}$$

where N is the number of the θ observations, θ (t_i) are the observations at a specific time, and θ(t_i,b) are the corresponding models' (Eqs. 1 to 3) predictions for the vector of optimized parameters, b (Ks and β). The inverse problem was solved using the fminsearch function in MATLAB. To evaluate the prediction quality, the root mean squared error (RMSE), the Nash–Sutcliffe efficiency (NSE), and the Pearson correlation (r) were calculated following Ritter and Muñoz-Carpena (2013).

3 Results and discussion

The water level measurements at the soil surface were collected by the Shafdan operators as part of the operational routine. Analysis of the water level data (ponding depth) indicates two operational modes of long and short cycles that were implemented at the investigated pond. The characteristics of the drying and wetting stages, as recorded by the operators, are summarized in Table 1. Note that the wetting stage is defined when water is present at the soil surface and the drying stage is defined when water is absent. Further note that the water level is measured at a single point (close to the pond inlet). Therefore, some variations, at the order of a few centimetres may exist due to microtopography and the distance between the inlet and the far parts of the pond. This may lead to some delays in water arrival or recession from the stations. Throughout the analysis described below, we define the winter period as the months between November and April and the summer period as the months between May and October, corresponding to the Mediterranean climate. Our monitoring systems are operated independently from the Shafdan facilities.

Table 1 Technical information of the recorded long and short wetting and drying cycles (these data were provided by the Shafdan operators).

3.1 Hydrological conditions

A representative set of θ time series measured in the SAT's vadose zone is presented to describe the variability in hydrological conditions measured throughout different seasons and operational modes (Fig. 2). The θ measurements were obtained during long (Fig. 2a, b, c) and short (Fig. 2d) cycles at three different depths. Note that the water content measurements presented in Fig. 2a and b were recorded at station 2 (Fig. 1) during summer (Fig. 2a) and winter (Fig. 2b). The water content variations under short cycles were measured at station 1 (Fig. 2c, d). Differences in the absolute values between the water content observations at different depths are mainly related to the vertical texture variability (Fig. S2). Under the long cycles, θ measurements were obtained throughout 19 cycles during summer (May–October) and 14 cycles during winter (November–April). Under short cycles, there were 12 cycles during summer and 25 during winter.

Every recorded wetting event prompted an intensive infiltration that was expressed by a rapid and almost instantaneous increase in water content at all depths and under the two operational modes (Fig. 2). Furthermore, the soil remained at similar level of saturation throughout each wetting stage. Similarly, once the drying process started, it occurred virtually simultaneously at all depths. There were noticeable differences in the drainage rates between summer and winter, where the soil dried faster in summer. To elaborate the drainage process, the drying stages were assembled and averaged at an hourly interval and separated into short (Fig. 3a) and long (Fig. 3b) cycles. Additionally, Eqs. (1), (2), and (3) (the hydrological model) were implemented to describe the water flow in the SAT's vadose zone under long and short cycles. The hydrological models were calibrated and validated against water content observations at 25 cm depth by adjusting the Ks and β parameters (Figs. S4 and S5; Table 2). Throughout the calibration, the Ks and β parameters attained different values for the long cycle periods and the short cycle periods (Table 2). There are differences in soil physical properties between stations 1 and 2, which explains the need for calibrating different parameter sets. In addition, the calibrated Ks values for both models were substantially lower than the measured Ks values (Fig. S2; Table 2). It has been shown that Ks measurements in the field are commonly lower than lab Ks measurements (Nimmo et al., 2009). This is related to a reduction in soil conductivity due to air trapping in the soil pores during the wetting process when water is applied at the land surface (Mizrahi et al., 2016; Nimmo et al., 2009).

Under the short cycles, the soil drainage process occurred mostly within the first 15 h of the drying stage (Fig. 3a). The soil drainage rate was slightly higher in summer than in winter. Under the short cycles, the model successfully followed the observed trends, where the validation period showed similar performances (Figs. S4, 3a). The model results for the short cycles confirm that the differences in drainage between summer and winter are mainly due to temperature changes that affect water viscosity (Lin et al., 2003). During the long cycle application, the drainage rates in summer showed a moderate θ decline compared to the observed θ under short cycles (Fig. 3b). This might be due to the differences in soil physical parameters between stations 1 and 2, as highlighted by the calibrated models' parameters (Table 2).

While the model under long cycles successfully followed the drainage trend during summer, the model showed poor performance during winter under long cycles. This is mainly due to the observable changes in θ measurements, which displayed a shift towards higher values from November 2019 (Figs. S5 and 2b). To explore the changes in the SAT's physical properties, an additional parameter set was calibrated against the winter data only during the long cycles (Fig. 3b; green line). Both the Ks and β parameters attained lower values, and the θ_s increased (Fig. S6; Table 2). Previous studies related the accumulation of organic matter in SAT to lower rates of organic matter decomposition during winter (Nadav et al., 2012, 2021a, b; Arye et al., 2011). The authors ruled out the occurrence of soil clogging and indicated that the accumulation of organic matter at the topsoil increased the degree of soil water repellency or soil hydrophobicity. This phenomenon often develops in sandy soils (commonly used in SAT) due to the low specific surface area of sand (∼ 0.0077 g m⁻²) compared to clay (∼ 900 g m⁻²; Doerr et al., 2000; Wallis and Horne, 1992). Thus, only a small amount of organic matter is required to coat the particles of the sand in order to develop soil water repellency (Wallis and Horne, 1992). Arye et al. (2011) showed that soil hydrophobicity is attributed to the reduction in the liquid surface tension and increase in the contact angle. These changes in soil properties are related to the reduction in the soil permeability (Nadav et al., 2012b). It appears that long wetting and drying cycles in SAT during winter can alter the physical soil properties, which eventually affect the infiltration capabilities.

Figure 2 Representative time series of θ measurements obtained at station 2 (Fig. 1) during (a) summer and (b) winter under long wetting and drying cycles at three different depths. The black line represents the surface water hydraulic head. The bottom plots display θ measurements obtained at station 1 under (c) long and (d) short wetting and drying cycles during the summer.

Table 2 Estimated parameters of the hydrological models.

Figure 3 The average and standard deviation values of measured θ at 25 cm depth throughout the drying stages at an hourly timescale. (a) Short cycles (station 1) and (b) long cycles (station 2) are shown. The blue and red circles represent the average θ values collected during winter (November–April) and summer (May–October), respectively. The statistics of measured θ under long cycles are based on 19 drying stages during summer and 14 drying stages during winter. For the short cycles, the statistics are based on 12 drying stages during summer and 24 drying stages during winter. The dashed black and solid grey lines represent the average values of simulated θ throughout the drying stage during winter (November–April) and summer (May–October), respectively. Note that, for the long cycle periods, an additional model was established based only on winter data (green line).

3.2 Seasonal differences in SAT redox (E_h) conditions under long wetting cycles

E_h conditions were monitored at four depths (Fig. 4a, b). The behaviour of the E_h dynamic following a wetting event was similar for the summer and the winter periods. At 25, 50, and 75 cm depths, a gradual decline in E_h started only after a time lag from the beginning of a wetting event (Fig. 4a, b). This time lag can be explained by the presence of dissolved oxygen (DO) in the percolating solution. Once the DO is depleted, suboxic and anoxic conditions begin to develop (Dutta et al., 2015; Ben Moshe et al., 2021, 2020). The E_h conditions at 25 cm were the most highly responsive to the wetting events, while the E_h conditions were the most negative at this depth (Fig. 4a, b). At 100 cm depth, only minor changes were observed, and in some cases (during winter), no changes were observed (Fig. 4a, b). According to the gaseous O₂ measurements (Fig. 4c, d), there were partially aerated conditions (unsaturated conditions) during some of the flooding events at 75 cm depth, while at 150 cm depth, unsaturated conditions prevailed continuously throughout the period of measurements. Therefore, the small changes in E_h at 100 cm depth could either have been the outcome of only minor biochemical activity or they could have been due to the sufficient oxygen supply. Gorski et al. (2019) suggested that, in coarse-grained soils, due to the high infiltration rates, the DO is delivered rapidly by the percolating water to the base of the saturated zone. Thus, the soil microbes do not have enough time to consume the DO. In the current study, it is reflected by positive E_h values at 100 cm depth (Fig. 4a, b). As the wetting cycle continues, the E_h conditions at 100 cm depth decrease, mainly during summer (Fig. 4a, b). This indicates that eventually the infiltrated water that reaches 100 cm depth contains lower concentrations of DO. Yet, the E_h conditions do not show a further decrease, which can be attributed to limited microbial activity. Previous studies showed that carbon availability at greater depths of the vadose zone is a dominant factor for this limitation (Grau-martínez et al., 2017, 2018; Brettar et al., 2002). Furthermore, the monitored E_h conditions illustrate that most of the activity in SAT systems occurs at the topsoil, as illustrated in previous studies (e.g. Quanrud et al., 1996, 2003; Sopilniak et al., 2017; Grinshpan et al., 2022). Note that the E_h recovery time, i.e. an increase towards positive values, was virtually instantaneous once the drying process initiated (Fig. 4a, b). The increase in E_h conditions occurred concurrently with the observed rapid increase of the gaseous oxygen (O₂) in the vadose zone (Fig. 4c, d).

A distinct difference in E_h conditions between winter and summer at 25 cm depth is expressed by a more negative range of values during summer (Fig. 4a, b). Since similar wastewater quality and hydraulic loads were fed to the pond during summer and winter (see the Supplement), the E_h conditions were mainly affected by the SAT's aeration state and seasonal temperature changes. In Fig. 5, the monthly E_h conditions at 25 cm depth are presented in the form of a box plot, together with the monthly average ambient temperature (dashed black line) and monthly average global radiation (grey line). The E_h conditions showed a wider range of values with the increase in temperature (Fig. 5). Between November and March, the E_h conditions mostly remained above zero or were slightly negative (Fig. 5). Once the temperature increased above 24 ^∘C between May and October, E_h conditions showed substantial fluctuation between negative and positive values (Fig. 5). Note that the average monthly ambient temperatures during November and May were similar in value, but during November, the E_h conditions mostly remained above zero (Fig. 5). This difference is connected to the daylight duration, as indicated by the global radiation, which is substantially greater in May than in November (Fig. 5; grey line). A typical characteristic of aquatic systems is the large fluctuations in dissolved oxygen (DO) concentrations due to intense photosynthesis and respiration (Stumm and Morgan, 1996). Goren et al. (2014) illustrated that following the effluent spreading in the infiltration ponds, photosynthesis enriches the water with DO. Furthermore, chemical analysis of porewater samples that were obtained at 0.5 m depth in the SAT vadose zone indicated the DO's substantial influence on the biochemical conditions of the percolating water (Goren et al., 2014). Thus, the photosynthesis process enriches the effluent with DO, which encourages further microbial activity (Goren et al., 2014; Hargreaves, 2006; Rodríguez-Escales et al., 2020). However, between July and September, there was a decrease in global radiation that did not affect the E_h variability (Fig. 5). Thus, it appears that, under long wetting stages, the seasonal temperature changes dominate the E_h conditions but show some trade-off with the global radiation.

Figure 4 Representative time series of redox potential (E_h) and gaseous O₂ measurements obtained at station 2 (Fig. 1) at depths of 25, 50, and 100 cm. (a, c) Summer and (b, d) winter results are shown. Note that the grey and white areas indicate wetting and drying periods, respectively.

Figure 5 Seasonal changes in E_h conditions as observed at 25 cm depth at station 2 under long (7 d) wetting periods. The dashed black line represents the monthly mean ambient temperature, and the solid grey line represents the monthly mean global radiation obtained from the Israeli Meteorological Service (IMS, 2022).

3.3 A comparison between SAT redox (E_h) conditions under long and short cycles

It has been shown that the operational mode affects the aeration conditions of the upper vadose zone, which, consequently, might alter the infiltration rates and the intensity of the biogeochemical processes (Goren et al., 2014). Throughout the measurement period, the intervals of wetting and drying stages were modified (Table 1). The wetting and drying stages were substantially shortened during September 2020, which enabled the examining of the differences in E_h conditions under short and long cycles. The variations in E_h conditions under long cycles during October 2019, December 2019, and January 2020 are presented in Fig. 6a and c. In addition, the changes in E_h conditions under short cycles during September–October 2020, December 2020, and January 2021 are presented in Fig. 6b and d.

As was shown above, under long cycles, the E_h conditions declined towards slightly negative values during winter and attained markedly negative values during the summer months at 25 and 50 cm depths, where only minor variations were observed at 100 cm depth (Fig. 6a, c). The E_h conditions, under short cycles during the summer, showed a decline towards slightly negative values for a brief time, compared to long cycles (Fig. 6a, b). During the winter, only minor variations in E_h conditions were recorded under short cycles (Fig. 6d). Note that, at 100 cm depth, there was almost no change in E_h conditions for either season under short cycles (Fig. 6b, d). As was illustrated previously (Fig. 4a, b), once the drying stage initiated, the recovery of E_h towards positive values (oxic conditions) was almost instant under both short and long cycles. It appears that the re-establishment of oxic conditions occurs independently of the length of the wetting stage.

The observations of E_h and gaseous O₂ under long and short cycles indicate a weak relationship between the wetting and the drying stages. These results call into question the advantage of implementing the wet $/$ dry ratio for optimizing SAT performance. It appears that the length of a wetting stage and a drying stage should be defined separately (Ben Moshe et al., 2020, 2021). Thus, the necessary further investigation is described below to examine the optimal lengths of the wetting and drying stages.

Figure 6 Representative time series of redox potential (E_h) measurements at 25, 50, and 100 cm depths at station 1 (Fig. 1). (a, c) Long wetting periods and (b, d) short wetting periods are shown. Note that the grey and white areas indicate wetting and drying periods, respectively.

3.4 The length of the wetting stage according to E_h measurements

During wetting stages, the main limiting factor of the biodegradation process is the availability of dissolved oxygen (Skopp et al., 1990; Cook and Knight, 2003). Once the soil pores are fully saturated, only the dissolved oxygen (DO) of the percolating water is available. However, studies have illustrated that the DO of the percolating water rapidly depletes (Dutta et al., 2015; Ben Moshe et al., 2020, 2021). Thus, the wetting stage should cease when the DO no longer has an effect on the degradation process, i.e. when suboxic and anoxic conditions begin to develop. E_h measurements can provide a good indication for the time when such conditions begin to be established. For this purpose, the hourly measured E_h conditions during the wetting stage of 33 recorded cycles (19 cycles during summer and 14 cycles during winter) were averaged for 25, 50, 75, and 100 cm depth (Fig. 7). The data were separated between winter (November–April) and summer (May–October), according to the trends presented in Fig. 5.

During summer, the decline in E_h conditions was steeper compared to the trends exhibited in winter (Fig. 7). Following 30 h of wetting, the E_h conditions during the summer dropped below 400 mV at 25 cm depth, which indicates the establishment of suboxic conditions (Fig. 7a). After 37 h, the E_h observations at 50 and 75 cm depths showed similar trends (Fig. 7a). Note that, after 45 h the E_h observations at 25 cm depth reach anoxic conditions. The differences between E_h measurements at topsoil and at deeper parts of the vadose zone were previously attributed to carbon availability, which decrease with depth (Brettar et al., 2002; Bohrerova et al., 2004). Furthermore, Brettar et al. (2002) provided a range of threshold E_h values for the denitrifying conditions, depending on the depth of the soil horizon. Similarly, despite the differences between the E_h measurements at 25 cm depth and at deeper depths (Fig. 7), all the E_h measurements potentially indicate similar activity. To prevent anoxic conditions in the SAT vadose zone during summer, the optimal length of a wetting stage, in terms of biodegradation, would be about 30 h (1.25 d). This is in agreement with Dutta et al. (2015), who suggested a relatively wide distribution of de-oxygenation times of between 0.36 and 1.5 d.

In winter, the E_h measurements decrease towards suboxic conditions, that occur after 51, 69, and 84 h of wetting at 25, 50, and 75 cm depths, respectively (Fig. 7b). The delay in E_h drop (compared to summer) is partly explained by the longer presence of gaseous O₂ in the SAT vadose zone (Figs. 4c, d and 8). It takes a longer time for the gaseous O₂ to deplete during winter (Fig. 8), probably as a result of lower oxygen demand due to slower microbial activity at low temperatures (Kirschbaum, 1995). An additional explanation for the longer presence of gaseous O₂ is due to a longer establishment of saturated conditions as a consequence of lower Ks values (see Sect. 3.1). However, once the gaseous O₂ is absent from the SAT vadose zone during winter (Fig. 8), no further decrease in E_h conditions occurs (Fig. 7b). Therefore, the trade-off between nutrient transport rates and reaction rates should be considered, where the supply of the percolated effluent's substrate might be faster than the SAT's degradation capability (Greskowiak et al., 2006). Determining the length of a wetting stage during winter, using solely the E_h and gaseous O₂ measurements, may be challenging. Nevertheless, the changes in soil physical parameters and the seasonal differences in E_h values imply low biodegradation rates during winter. From a practical perspective, the observed E_h values at 25 cm depth during winter are compared with the E_h measurements during summer. Using a two-sample t test shows that the differences in the first 28 h of a wetting cycle are insignificant. It suggests that the length of a wetting stage during winter should be no more than 30 h, as in the summer period. Nevertheless, given the evidences that are reflected in the measurements of the current study and conclusions from previous studies, the implementation of SAT during winter is questionable. It might be necessary to store the effluents in reservoirs during winter and apply the SAT during summer.

According to the E_h sensor at 100 cm, suboxic conditions prevail for most times during winter and a moderate decrease (compared to the other sensors) in E_h conditions is noticeable during summer (Fig. 7). The gaseous O₂ observations at 75 cm indicate depletion of gaseous O₂ during summer and a very low O₂ concentration (2 %) during winter (Fig. 8c). At 150 cm, the gaseous O₂ observations suggest that, for most of the recorded wetting cycles, unsaturated conditions prevailed (Fig. 8d). It is not clear if the E_h conditions at 100 cm are determined by a continuous supply of oxygen, given the unsaturated conditions. This might affect the decision concerning the optimal length of a wetting cycle. However, previous studies (Miller et al., 2006; Fox et al., 2005; Lin et al., 2008; Goren et al., 2014; Sopilniak et al., 2018; Essandoh et al., 2013; Quanrud et al., 1996, 2003) and a recent study by our group (Grinshpan et al., 2022) suggest that most biodegradation activity occurs at topsoil, and there is a steep reduction in removal capabilities with depth. Thus, the potential contribution of the deeper vadose zone to the SAT treatment should be further investigated.

Figure 7 Plot of hourly means of E_h measurements obtained at station 2, from the beginning of the wetting stages, at three depths. The statistics are based on 19 cycles during summer (May–October) and 14 cycles during winter (November–April) that were observed in (a) winter and (b) summer. Note that the data are separated between summer and winter according to the trends presented in Fig. 5. The green, grey, and yellow areas represent the oxic (400 mV <), suboxic (400 mV > E_h > 100 mV), and anoxic E_h < 100 mV) conditions, respectively.

Figure 8 Plot of hourly means of gaseous O₂ measurements from the beginning of the wetting stages at (a) 25, (b) 50, (c) 75, and (d) 150 cm depths. The statistics are based on 14 cycles during winter (November–April) and 19 cycles during summer (May–October) that were observed at multiple depths in station 2. Note that the data are separated between summer and winter according to the trends presented in Fig. 5.

3.5 The length of the drying stage using E_h and gaseous O₂ measurements

The drying stage in SAT systems is implemented to restore the infiltration, biological, and chemical capabilities of the pond, mainly by aerating the soil (Sharma and Kennedy, 2017). A drying stage is defined as the stage in which there is no water at the soil surface, and the observed volumetric water contents begin to decrease (Fig. 2). To determine the optimal time for a drying stage, both the averaged values of the E_h (Fig. 9a, b) and gaseous O₂ (Fig. 9c, d) observations during the drying stages are presented. Note that the E_h measurements were conducted at 25, 50, 75, and 100 cm depths, and the O₂ observations were obtained at 25, 50, 75, and 150 cm depths (Fig. 9). The E_h and O₂ data during the drying stages were separated between winter (November–April) and summer (May–October), according to the trends presented in Fig. 5.

Slower aeration rates and larger variability were observed during winter compared to summer (Fig. 9a, b). Furthermore, the E_h recovery is clearly dominated by the rates of the gaseous O₂ intrusion to the soil (Fig. 9c, d). During winter, the aeration rates are different at each depth, but in general, the gaseous O₂ concentrations increase moderately with time (Fig. 9a). The recovery of the E_h conditions display comparable trends (Fig. 9c). Although the O₂ sensor at 50 cm depth suggest a faster O₂ intrusion rate compared to the other sensors (Fig. 9a), it is not expressed in the E_h recovery time (Fig. 9c). During summer all the observed O₂ curves show steep recovery that ceases after about 20 h (Fig. 9b). According to the O₂ observations, the gaseous O₂ intrusion process at 50 and 75 cm depths starts at the very beginning of the drying stage, while at 25 and 100 cm depths there is a 2 h delay (Fig. 9b). Nevertheless, the E_h conditions recovery is comparable at all depths (Fig. 9d).

The E_h observations illustrate that the re-establishment of oxic conditions in the SAT vadose zone is similar during winter and summer (Fig. 9c, d). Once the drying stage started, it required about 36 h for the E_h conditions to reach values in the range of oxic conditions, regardless of the initial value of the E_h conditions (Fig. 9c, d). However, a recent study has indicated that using solely the gaseous O₂ concentration to quantify soil aeration status might be inaccurate for some conditions (Ben-Noah et al., 2021). Instead, the authors suggested using the Damköhler number, which is the ratio between characteristic diffusion (i.e. oxygen supply) and soil respiration times. A small Damköhler number would indicate that the O₂ consumption rate is lower than the O₂ supply. Thus, although during summer the oxygen supply rate is relatively high compared to winter (Fig. 9a, b), the oxygen consumption rate for soil respiration is expected to be substantially higher during summer than during winter (Kirschbaum, 1995), as it is reflected by the E_h observations (Figs. 5 and 7). Therefore, in parallel to proper treatment of aeration time, optimal wetting time should be further treated.

Figure 9 Plot of hourly means of E_h and gaseous O₂ measurements from the beginning of the drying stages. The statistics are based on (a, c) 14 cycles during winter (November–April) and (b, d) 19 cycles during summer (May–October). Note that the gaseous O₂ (a, b) measurements were obtained at 25, 50, 75, and 150 cm depths, and the E_h conditions (c, d) were observed at 25, 50, 75, and 100 cm depths. All observations were conducted in station 2. The data are separated between summer and winter according to the trends presented in Fig. 5. The green, grey, and yellow areas represent the oxic (400 mV <), suboxic (400 mV > E_h > 100 mV), and anoxic E_h < 100 mV) conditions, respectively.

4 Summary and conclusions

Continuous monitoring of E_h, θ, T, and gaseous O₂ in the vadose zones of SAT infiltration ponds was carried out for about 600 d. SAT operation was subjected to long and short wetting and drying cycles and seasonal changes. The datasets enabled the examination of factors that control the hydrological and geochemical conditions in SAT. Calibrated and validated hydrological models were applied to investigate the water flow dynamics in the SAT vadose zone. E_h and gaseous O₂ observations were averaged to determine seasonal changes and deduce the optimal length of wetting and drying stages.

The examination of the measured E_h conditions illustrated a noticeable decline to markedly negative values (−450 mV >) during summer and to values between 0 and −50 mV during winter. These E_h conditions were established following 30 h of wetting, and no considerable changes in E_h conditions were noticeable until the wetting stages ceased. A monthly statistic of the E_h measurements illustrated the relationship between the size of the E_h amplitude and the seasonal T changes. Furthermore, it is speculated that the limited decrease in E_h conditions during winter was due to lower microbial activity. An additional support for this claim is the reduction in infiltration capabilities following long wetting and drying cycles during winter, as was indicated by the hydrological models. To define the optimal length of a wetting stage, the E_h data were averaged and separated between the winter and summer periods. During the summer period, the optimal time length of a wetting stage, according to E_h observations, is about 30 h. Determining the optimal length of a wetting stage during winter is challenging. Practically, there are no significant differences between E_h conditions during winter and summer for the first 28 h of wetting. Thus, the length of a wetting stage during winter should be 30 h, as in the summer period.

The length of a drying stage, following the E_h observations, should be about 36 h, regardless of the initial values of the E_h conditions and the season. Note that, during summer, a longer drying time is required for the E_h conditions to attain suboxic conditions. Nevertheless, the gaseous O₂ observations indicated faster aeration rates during summer, which compensate for the very negative E_h conditions and allow fast recovery. Following our analysis, we suggest that, under the tested conditions, the length of a drying stage should be 36 h for winter and summer. Shorter drying stages would affect SAT efficiency, while applying longer drying stages would reduce the total hydraulic loads that can be fed to the infiltration basin.

Implementing in situ E_h, θ, T, and gaseous O₂ sensors provided continuous high-resolution observations. These datasets revealed the hydrological and biochemical dynamics in SAT as imposed by seasonal and operational changes. Analysis of the E_h and O₂ measurements enabled the identification of the optimal time lengths of wetting and drying stages. The results indicated that there is no direct relationship between the length of a wetting stage and a drying stage. Thus, the operational use of a wetting/drying ratio in SAT management should be reconsidered.