Evaluation of hydrodynamic characteristics of porous media from one-step outflow experiments using RETC code

The knowledge of two main hydraulic properties, water retention curve θ(h) and unsaturated hydraulic conductivity Κ(θ), is necessary in studying the movement of water and soluble salts in unsaturated soils, as well as in their simulation models. The determination of soil hydraulic properties can be achieved either with field-based or laboratory methods. Each method is characterized by specific properties with their limitations and inherent assumptions (Bordoni et al. 2017). However, in laboratory methods, while the water retention curve θ(h) can be determined relatively quickly and easily, the determination of hydraulic conductivity as a function of soil water content, K(θ), or pressure head, Κ(h), is a difficult and time-consuming process. For this reason, several statistical models of pore size distribution have been developed for the indirect prediction of K(θ) using water retention curve and saturated hydraulic conductivity K_s (Childs & Collis George 1950; Burdine 1953; Mualem 1976). The introduction of analytical expressions of the water retention curve θ(h) in combination with the abovementioned models led to closed-form analytical predictive models of hydraulic properties (Brooks & Corey 1964; van Genuchten 1980). Typically, the absence of experimental data, especially K(θ) data, has led to the widespread use of closed-form analytical predictive models in the prediction of hydraulic properties (Ghazanfari et al. 2016; Arrey et al. 2018). This suggests the necessity of assessing the accuracy of the hydraulic property predictions since unsaturated hydraulic conductivity plays an integral role in determining the accuracy of any numerical solution to water flow and contamination problems (Yates et al. 1992). In some cases, it has already been proved that the calculated values of K(θ) deviate significantly from measured values (Talsma 1985; Yates et al. 1992; Valiantzas et al. 2007; Londra 2010).

The one-step outflow method (Doering 1965) is one of the most widely used laboratory methods for determining the soil–water diffusivity relationship D(θ). The method can be easily adopted for routine laboratory work, e.g. is easily applied in the same soil sample and apparatus that are used for the determination of the soil–water diffusivity relationship D(θ), and the soil–water retention curve. The K(θ) relationship can then be calculated using D(θ) and water retention data (Childs & Collis George 1950).

Many researchers have proposed analytical methods for calculating D(θ) from the one-step outflow data without the required assumptions of any mathematical form for the hydraulic properties, which is an advantage (e.g. Passioura 1976; Valiantzas 1989; Valiantzas et al. 2007; Londra & Valiantzas 2011). However, in these cases, independently measured θ(h) data to determine unsaturated hydraulic conductivity K(θ) is required. Valiantzas et al. (2007) proposed a direct method based on a simple curve-fitting procedure applied to the experimental outflow data for a simple power and extended power form function which leads to direct calculation of the soil water diffusivity function from explicit formulae.

Alternatively, in many cases, parameter estimation methods together with outflow experiments for determining hydraulic properties have been used (e.g. Parker et al. 1985; van Dam et al. 1992; Bitterlich et al. 2004). However, an assumption of particular mathematical forms for the hydraulic properties is required.

The main disadvantage of the one-step outflow method is the weakness to predict K(θ) values near saturation, due to the fact that the outflow method cannot be applied at the first stage of outflow, where the flow rate is essentially controlled by the porous plate resistance (Passioura 1976). This disadvantage may be overcome using RETC code (van Genuchten et al. 1991) as has been studied recently in a sand sample by Bourazanis et al. (2016). More specifically, the comparison between K(θ) values near saturation measured by steady-state laboratory method and predictions obtained by RETC code using simultaneously experimental θ(h) values and K(θ) from one-step outflow method showed very good agreement.

The RETC code enables the calculation of the conceptual models parameters either by fitting only the water retention data for predicting the hydraulic conductivity function or by fitting simultaneously both water retention and hydraulic conductivity experimental data, assuming various analytical forms of θ(h) and K(θ), by including the effect of pore connectivity, i.e. parameter p of the Mualem model (Mualem 1976). Many researchers have used RETC code with different scenarios and number of model-fitting parameters with various results (van Genuchten & Leij 1992; Yates et al. 1992; Schaap & van Genuchten 2006; Kargas & Londra 2015; Bourazanis et al. 2016). However, in cases where the water retention and hydraulic conductivity experimental data are simultaneously used, it is rarely reported whether the experimental data were obtained from the same soil sample using the same apparatus. In addition, among the different methods of determining K(θ) relationship, there is a large difference in the range of water content and pressure-head measurements (Stolte et al. 1994). Also, Siltecho et al. (2015) demonstrated that the van Genuchten unsaturated soil parameters were significantly different according to the measurement methods employed.

These factors may play a significant role in the results of various prediction scenarios and the necessary number of fitting parameters.

From the research published so far, it seems that the prediction of soil hydraulic properties has not been considered when experimental θ(h) and K(θ) data are used simultaneously as input data in the RETC program when the K(θ) data have been obtained from the one-step outflow method. More specifically, when the K(θ) data are obtained by using the simplified equation of Valiantzas et al. (2007).

The main objectives of the present study were to experimentally determine, on the same soil sample and in the same apparatus, θ(h) and K(θ) relationships in six porous media using the one-step outflow method and then compare the experimental values with those predicted by assuming the Mualem-van Genuchten model (Mualem 1976; van Genuchten 1980) using RETC code (van Genuchten et al. 1991) for two prediction scenarios with three model fitting parameters (α, n, θ_r) each. In the case of the second scenario applied, the Mualem model parameter p as an additional fitting parameter was also examined.

Experiments to determine the water retention curves followed by the one-step outflow procedure to determine D(θ) and then K(θ) were performed in the laboratory for six disturbed porous media with different soil textures: a sand (0.2 mm < d < 0.5 mm), a sandy loam (13.2% clay, 8% silt, 78.8% sand), a loam (20% clay, 38% silt, 42% sand), a clay loam (28% clay, 37% silt, 35% sand), a silty clay loam (36.5% clay, 52% silt, 11.5% sand) and a clay soil (47% clay, 36% silt, 17% sand). Note that the clay soil aggregate was used with aggregation fraction 0.5–1.2 mm.

Water retention curve measurements and one-step outflow experiments were performed in the laboratory using: (i) a Haines-type assembly (Haines 1930) for the sand and sandy loam soil; and (ii) Richards' pressure cell chambers (Kargas & Londra 2015) for the loam, clay loam, silty clay loam and clay soil. Initially, the water retention curves θ(h) were measured, followed by the one-step outflow experiment.

Disturbed soil samples of the sand and sandy loam soil, 2.5 cm high and 9.6 cm diameter, were placed on a tension plate apparatus in a Haines-type assembly, and samples of the loam, clay loam, silty clay loam and clay soil, 3 cm high and 7 cm diameter, were placed in a Richards' pressure cell chambers.

Soil samples were allowed to wet from the bottom of the tension plate gradually until saturation. After that, the samples were subjected to a drying-wetting cycle and the primary drying water-retention-curve data were obtained by applying: (i) negative pressure steps through the saturated tension plate, in the case of Haines-type assembly; and (ii) gas pressure steps to the top of the soil sample, in the case of Richards' cell, weighting the water lost at various pressure steps.

At the end of the water-retention-curve measurement, on the same samples, in the same apparatus, saturation of the samples was performed, followed by the one-step outflow procedure.

A large negative or positive pressure step h_f (L) was suddenly applied at the bottom of the sample (Haines-type assembly) or at the top of the sample (Richards' cell), respectively, equal to the final pressure step used in the determination of the water retention curve (h_f = −90 cm for sand, h_f = −134 cm for sandy loam soil, h_f = +700 cm for loam, clay loam, silty clay loam and clay soil) and the cumulative outflow volumes V_i (L³) were recorded with time t_i (T) (i = 1, 2, 3 ... Ν). The time steps used were 5 seconds at the beginning of the experiment and then they were adjusted according to the outflow rate, depending on the soil type, reaching gradually 1,800 seconds at the end of the experiment. The total times were 6,450, 6,696, 12,098, 16,052, 18,372 and 24,035 seconds for sand, sandy loam, clay loam, loam, clay and silty clay loam soil, respectively.

In Figure 1, an example of schematic illustration of the experimental function for silty clay loam soil, and the three-parameter power fitting function (Equation (2)) is presented. The three stages of outflow, according to Passioura (1976), are clearly identifiable. In the first stage, corresponding to the initial part of the curve before the outflow becomes linear with ⁠, it is the plate impedance only that determines the outflow rate. In the second stage, during which the outflow varies linearly with ⁠, the effect of plate impedance has not yet become negligible. The third stage corresponds to the portion of the curve where cumulative outflow ceases to be linear with and the effect of plate impedance is minimal.

The slope (L⁻¹) was calculated from the experimental water retention curve.

The saturated hydraulic conductivity, Κ_s (LT⁻¹) was independently determined by the constant-head method (Klute & Dirksen 1986).

The model-fitting parameters described above were evaluated by RETC program from measured water retention and hydraulic conductivity data. The unknown parameters of the Mualem-van Genuchten (M-vG) model in the parameter optimization process to fit the water retention and unsaturated hydraulic conductivity functions were θ_r, α and n using two scenarios. In the first scenario, the experimental data of water retention curve and saturated hydraulic conductivity (K_s) were used as input data. In the second scenario, the experimental values of θ(h), Κ_s and K(θ) calculated from D(θ) values using experimental outflow data were used as input data and two types of conductivity model were examined for regression analysis. Type I: ‘Conductivity versus Water Content’ and Type II: ‘Logarithmically Transformed Conductivity versus Water Content’. The value of the parameter p was either taken as constant and equal to 0.5, a value widely used (Mualem 1976), or as a fitting parameter. Evaluation of the results was performed by comparing experimental and predicted values using root mean squared errors (RMSE).

In any case studied, the M-vG model fitting parameters α, n, θ_r, as well as the experimental values θ_s, Κ_s for all porous media examined, are given in Table 1. It is worthy of note that high predicted θ_r values of the fine textured soils examined can be attributed to the fact that they come from mathematical fitting which is affected by the low pressure step (h_f = +700 cm) applied. In this pressure, a sufficiently large proportion might not be drained from the total pore space. This indicates the need to apply a higher pressure step for measuring the hydraulic properties and evaluating the fitting parameters. Also, the fact that θ_r values of clay soil are lower than those of clay loam and silty clay loam should be attributed to the nature of the sample. As presented in the materials and method section, clay sample is not a typical soil like others but includes only an aggregated fraction of 0.5–1.2 mm soil particles.

As shown in Figure 2, it is apparent that there is a very good agreement of the results between experimental and predicted values of θ(h) for the first scenario, indicating that the corresponding soil hydraulic parameters listed in Table 1 provide an adequate description of θ(h) relationship. This is demonstrated by taking into consideration the acceptably small values of RMSE listed in Table 2.

On the other hand, the soil hydraulic parameters of the first prediction scenario, which is the most common used, using as input data only the experimental θ(h) and Κ_s data, do not adequately describe the K(θ) relationship (Figure 3). The K(θ) curves show a variable and unpredictable behavior, suggesting a non-reliable estimation based only on retention data θ(h). This is also shown by the large values of RMSE (Table 2).

In Figure 3 a comparison between Κ(θ) values from experimental one-step outflow data and predictions obtained by RETC for the two prediction scenarios are presented. The values of Κ(θ) were obtained by Equation (4) taking into account experimental one-step outflow data and water retention curve.

As shown in Figure 3, in the case of the second prediction scenario, using as input data simultaneously the experimental θ(h), K(θ) and K_s values for Type I of conductivity regression analysis, the prediction of K(θ) is satisfactorily improved compared to the first prediction scenario except in the case of clay soil, while for Type II of conductivity regression analysis, a very satisfactory prediction of K(θ) is obtained for all porous media studied. More specifically, better improvement in the prediction of Κ(θ) was observed in the case of Type II (logarithmically transformed conductivity) compared with both Type I (untransformed conductivity values) and the first scenario (Figure 3). This is also demonstrated in the RMSE values presented in Table 2. There is a decrease in the RMSE values indicating the abovementioned improvement. This effect is more pronounced for Type II.

Similar results were presented by Yates et al. (1992) which recommended that the logarithmically transformed values should be used since they better describe the conductivity across the entire range of observed values.

However, the cost of this improvement is a poorer characterization of θ(h) relationship. Although Type II of conductivity regression analysis accurately describe K(θ), it reduces the predictive accuracy of θ(h) compared with the first prediction scenario for all porous media studied. This can also be seen from the results in Table 2 where, for the water content, there is an increase in RMSE values.

In order to remove this weakness, we investigated the case of increasing the model fitting parameters. We studied the Mualem model parameter p as an additional fitting parameter using RETC for Type II. As shown in Table 3, the fitting parameter p, in all porous media studied, had values with remarkable deviation from the conventional value of 0.5 ranging from −0.591 to 2.622 resulting from RETC code analysis. Similar behavior has also been noticed by other researchers (Leij et al. 1992; Yates et al. 1992; Kargas & Londra 2015).

Taking into consideration the RMSE values shown in Table 3, it appears that Type II method with four fitting parameters compared with the corresponding one with three fitting parameters provides an improvement in the description of the water retention curve for all porous media studied (except clay soil). On the other hand, increasing of fitting parameters provides an improvement (lower RMSE values) of hydraulic conductivity or leads to similar results to those obtained by Type II with three fitting parameters.

Similar results were obtained by Yates et al. (1992) when they applied two prediction scenarios using as input data: (i) the experimental values of θ(h) and Κ_s with four fitting parameters, and (ii) simultaneously the experimental values of θ(h) and Κ(θ) for the two types of conductivity regression analysis (logarithmically transformed conductivity–untransformed conductivity) with five or six fitting parameters. More specifically, the simultaneous method with five or six parameters was found to be better, but there were only minor differences between two types of simultaneous methods.

Additionally, from the results we can reasonably assume that the use of RETC using simultaneously experimental θ(h) and K(θ) data from outflow measurements using the Valiantzas et al. (2007) equation for Type II conductivity regression analysis, may improve the description of the K(θ) near saturation. In this way, the weakness of one-step outflow method predicting K(θ) at saturation may be overcome. This approach seems to be especially suited for large-scale studies that require realistic simulations in the wet region, e.g. during infiltration into soils.

In six porous media with different soil texture, two basic hydraulic properties, water retention curve and hydraulic conductivity were determined on the same sample using an easy and fast methodology. The experimental K(θ) values obtained by one-step outflow data were compared with those predicted by the Mualem-van Genuchten model using RETC code for two prediction scenarios with three fitting parameters a, n, θ_r and two types of conductivity regression analysis (Type I, untransformed conductivity; Type II, logarithmically transformed conductivity). From the results of this group of soils, it appears that the first scenario using only the experimental θ(h) and K_s values supplies soil hydraulic parameters that provide an adequate description of θ(h) relationships compared to the other methods studied. However, this scenario is unable to adequately predict K(θ) relationships. In the case where RETC with input data the experimental θ(h), K(θ) and K_s values for Type II was used, the obtained K(θ) predictions were in very good agreement with experimental values. However, this method causes a decreasing agreement between experimental and predicted θ(h) relationship compared with the method of the first prediction scenario. The best method for predicting both the θ(h) and K(θ) relationships, for this group of soils, is to use Type II regression analysis with p as an additional fitting parameter compared to the other methods examined. Overall, the combination of RETC with the one-step outflow method is a powerful tool which may be used successfully as a routine procedure for the determination of K(θ) even in the case of K(θ) values near saturation.