Closed‐form hydraulic conductivity equations for multimodal unsaturated soil hydraulic properties

Closed‐form expressions of the hydraulic conductivity function for linearly superposed subretention (multimodal) functions were derived for arbitrary sets of the Brooks and Corey (BC), van Genuchten (VG), and Kosugi (KO) water retention models. The generalized Mualem hydraulic conductivity model was evaluated using the mathematical approach of Priesack and Durner. Three types of modification to the multimodel were also proposed. Firstly, the derived conductivity equations can be simplified when the submodel parameters, hbi${h_{{{\rm{b}}_i}}}$ for the BC model, αi−1${{\alpha}}_i^{ - 1}$ for the VG model, and hmi${h_{{{\rm{m}}_i}}}$ for the KO model have the same (common) value (denoted as CH). Secondly, as in the case of the modified single VG and KO models, a hypothetical air‐entry head near saturation can be introduced for the multimodal VG and KO models to prevent unrealistic reductions in the hydraulic conductivity near saturation when the VG n parameter approaches its lower limit of n = 1. Furthermore, the multimodal hydraulic conductivity functions become a simple sum of conductivity subfunctions when the exponent r is unity (such as for Burdine's model), which leads to independent tortuosity effects for each submodel. The models are illustrated for two soils: a highly aggregated Kumamoto Andisol and a relatively unimodal dune sand. The dual‐(BC, VG, KO) and the VG1BC2 models equally represented the water retention data of the Andisol, with similar hydraulic conductivity curves. The dual‐BC‐CH, dual‐VG‐CH, and VG1BC2‐CH models fitted the water retention data of the dune sand well, with the hydraulic conductivity curves of the dual‐porosity model being similar to those of the Fayer and Simmons (FS) model.


INTRODUCTION
Unimodal functions of the soil hydraulic properties, such as the standard Brooks and Corey (1964) and van Genuchten (1980) models, are not appropriate for simulating unsaturated flow in macroporous soils and fractured rock. Durner (1994) proposed a multimodal water retention function (WRF) for these media by linearly superposing the van Genuchten (VG) model as subretention functions. Although Seki (2007) similarly defined a multi-Kosugi model (Kosugi, 1996) to fit various soil water retention data, analytical expressions of the hydraulic conductivity function (HCF) for the Mualem (1976)-based multi-WRFs were not derived. Priesack and Durner (2006) showed a procedure to evaluate the  Mualem (1976) model for linearly superposed subretention functions and derived a specific expression of the HCF using the multimodal VG model. Their closedform multi-VG model has been widely used to simulate water flow in macroporous and aggregated soils (Dimitrov et al., 2014;Lipovetsky et al., 2020;Schelle et al., 2010;Watanabe & Osada, 2016). Assouline and Or (2013) showed that the impact of the subsystem on the multimodal HCF can be significant, even when the subretention function has a relatively minor effect on the multimodal WRF. It is hence important to carefully select appropriate submodels for the multimodal formulations.
Alternatively, several WRFs accounting for capillary and adsorptive retention have been proposed. For example, Fayer and Simmons (1995) (the FS model) expressed the main WRF of the capillary region with the Brooks and Corey (1964) model (the BC model) or the VG model, while expressing the adsorption region with the Campbell and Shiozawa (1992) equation, where the water content decreases exponentially to 0 at some large value (h m ) of the pressure head, h, instead of using the residual water content, θ r . Groenevelt and Grant (2004) proposed using a relative humidity of 1% for anchoring the value of ℎ m (8.2 × 10 6 cm). Fayer and Simmons (1995) evaluated the integral of Mualem's equation in their model using power series expansions, leading to a relatively complicated closed-form expression for the HCF with many subfunctions. Peters (2013) further proposed a conductivity model accounting for the contributions of capillary, film, and vapor conductivity from saturation to zero soil water content at oven-dry conditions. Linear superposition of the WRF was used for the capillary subfunction given by the VG or KO models, and a film subfunction based on the BC model. Whereas Mualem's estimation of the HCF was used only in the capillary range, one additional parameter independently determined from the WRF was used to quantify the contribution of film conductivity. Since the method of deriving the HCF of linearly superposed WRFs (Priesack & Durner, 2006) is general and can be applied to other subfunctions such as the BC and KO models, the closed-form hydraulic functions based on Mualem's conductivity model over the whole range of capillary and the film flow would be very useful for many practical flow simulations.
In this note, closed-form equations of the HCF based on the generalized Mualem's hydraulic conductivity model are derived for linearly superposed subretention functions using the mathematical approach of Priesack and Durner (2006). The multimodels consist of an arbitrary number and combination of the BC, VG, and KO submodels. Three types of modified equations are also presented for giving flexibility in the application of the multimodel. Application of the dual submodels is demonstrated with literature retention datasets for aggregated volcanic and a coarse-textured medium.

Core Ideas
• Mualem-based hydraulic conductivity functions for multimodal retention models are derived. • The approach uses the Brooks-Corey, van Genuchten, and Kosugi water retention models. • The multimodal equations describe the hydraulic properties from saturated to dry conditions. • The closed-form equations are attractive for flow simulations because of their simple form.

HYDRAULIC CONDUCTIVITY FUNCTIONS
Effective saturation S as a function of the pressure head h is defined by where θ is the volumetric water content, and θ r and θ s are the residual and saturated water contents, respectively. Durner (1994) defined S(h) of a multimodal soil system as the sum of subsystem effective saturations S i (h): where k is the number of subsystems, and w i are weighting factors with 0 < w i < 1 and Σw i = 1. For notational convenience, we use the pressure head, h, to be positive for unsaturated conditions, thus effectively considering h to be an equivalent suction. Many formulations can be used for S(h) in Equation 1, and/or for the subsystem saturations S i (h) in Equation 2, such as the BC, VG, and KO functions. A similar linear combination, however, is not possible with most statistical poresize distribution models for predicting the hydraulic conductivity from available water retention data, such as those by Burdine (1953) and Mualem (1976). A generalized form of predictive conductivity models is given by (Assouline & Or, 2013;Kosugi, 1999;Kosugi et al., 2002;Priesack & Durner, 2006) where K is the unsaturated hydraulic conductivity, K s is the saturated hydraulic conductivity, K r is the relative hydraulic conductivity, and where p, q, and r are constants (with q > 0 and r > 0). In particular, p = 0.5, q = 1, and r = 2 for Mualem's model, and p = 2, q = 2, and r = 1 for Burdine's model. Priesack and Durner (2006) showed that the integral in Equation 3 for the multiple retention function can be given as a summation of the integral for the subretention functions using the chain rule: The overall S(h) of Equation 2 can be any arbitrary combination of subretention functions given by BC, VG, and/or KO models. The BC subretention function is expressed as where ℎ b > 0 is the air entry head and λ > 0. The VG subretention function is expressed as where α > 0, 0 < m i < 1, and n i > q for the parameter q > 0 as used in the generalized conductivity Equation 3. When VG-Mualem's hydraulic conductivity model is used, q = 1 and m i = 1 − 1/n i , whereas in the VG-Burdine model q = 2, and hence m i = 1 − 2/n i (van Genuchten, 1980). The KO subretention function is given by where ℎ m > 0, σ > 0, and Q is the complementary cumulative normal distribution function defined by , the HCF can be described in terms of subretention functions S i (h) as The closed-form HCF can be obtained by substituting the subretention function (ℎ) of Equations 5-7 into Equation 9 and evaluating the integral using the same technique as used for the subhydraulic conductivity function. Table 1 presents (h) and for three sub-retention functions (ℎ). We note that Priesack and Durner (2006) missed the exponent of α in (h) and in their generalized HCF for the multi-VG (Equations 5 and 6 in Priesack & Durner, 2006), but this does not affect Mualem's model ( = 1).
For a single subretention function, k = 1, in which case Equation 9 reduces to and substituting 1 (ℎ) and 1 (Table 1) in Equation 12 leads then immediately to the familiar forms of Mualem-type hydraulic functions for the BC, VG, and KO models. The combined model is denoted here by subscripting the number of the subfunction. For example, VG 1 BC 2 model denotes VG-type for S 1 (h) and BC-type for S 2 (h). Combinations of the same subfunctions (e.g., BC 1 BC 2 BC 3 . . . ) are referred to as the multimodels (e.g., multi-BC). The multi-VG model is the same as used in Durner (1994) and Priesack and Durner (2006). Multimodels consisting of only two similar subfunctions are referred to as dual-models, such as the dual-BC for BC 1 BC 2 and the dual-VG for VG 1 VG 2 . As an example, the VG 1 BC 2 model is given by   Table 1 become − for all the subsystems. These common terms in and are then cancelled in Equation 9. We use CH (common head) for this assumption, such as the multi-BC-CH model for which = ℎ b = ℎ b , and the VG 1 BC 2 -CH model for which = α −1 1 = ℎ b 2 . In case of the multi-VG-CH model with = α −1 , the denominator becomes ∑ =1 = 1 since is − regardless of and the − terms in the denominator and numerator are cancelled. Hence r (ℎ) can then be simply written as where (ℎ) = [1 + (αℎ) ] − . As will be discussed in Section 3 for Hamaoka dune sand, this assumption may be useful for soils with a unimodal pore-size distribution to reduce the number of parameters. Secondly, Equation 9 can be used also for the modified forms of the multi-VG and multi-KO models. Since the VG model shows unrealistically large reductions in the hydraulic conductivity near saturation when n approaches its lower limit of 1.0, Vogel et al. (2000) introduced a hypothetical air-entry head h b near saturation (such as h b = 1 cm) to stabilize flow simulations. Kosugi (1994) also proposed a modified form in the same way. Modified multi-VG and -KO models result by replacing Equation 2 with the following WRF: where h b is the air-entry head, and θ m a fictitious water content at h = 0, giving θ s at h = h b [i.e., S(h b ) = 1]. Using the same procedure as for the multi-VG and KO models, K r (h) can be calculated as where S i (h) and A i (h) are the same as those for the regular VG and KO models in Table 1. Notice that Equations 16 and 17 are identical to those given by Vogel at al. (2000) and Kosugi (1994) in case of k = 1. Finally, as Priesack and Durner (2006) pointed out, for r = 1, which includes Burdine's model (p = 2, q = 2, r = 1), the multimodal unsaturated conductivity function of Equation 9 can be written as a sum of conductivity subfunctions: where γ = ∑

=1
. Since the superposition of conductivity subfunctions becomes now linear, the parameter in this case can be specified independenty for each submodel: which shows that one can assign values to give independent tortuosity effects on the hydraulic conductivity.

APPLICATIONS
The derived models are illustrated here for published data of a Kumamoto Andisol (Sakai & Toride, 2007) and a Hamaoka dune sand (Miyazaki, 1976;Sakai et al., 2009) in Japan. The Kumamoto Andisol was selected as an example of a highly aggregated soil, whereas Hamaoka dune sand exemplified a soil having unimodal pore structure and adaptive retention properties. The parameters for the selected WRF models were estimated from measured (h i , θ i ) data using least-square T A B L E 2 Estimated water retention function parameters of Kumamoto Andisol (Figure 1) and Hamaoka dune sand (Figure 2) soils, and the coefficients of determination (R 2 ) of the regression Note. θ s , saturated water content; w 1 , weighting factor; ℎ b and λ i are parameters for the BC subretention function; α i and n i are parameters of the VG subretention function; optimization. As discussed further below, since the multimodels can describe water retention data in the very dry range using the second or higher submodels, we fixed θ r at 0, in which case θ(h) simplifies to θ s S(h), leading to a reduction in the number of parameters. The selected models are then the dual-(BC, VG, KO) models and the VG 1 BC 2 model for Kumamoto Andisols, and the dual-BC-CH, dual-VG-CH, VG 1 BC 2 -CH, and FS (Fayer & Simmons, 1995) models for Hamaoka dune sand. The estimated parameters and coefficients of determination of the least-square regression are shown in Table 2.
The HCFs were drawn using the determined WRF parameters for Mualem's model (q = 1, r = 2), with p = 0.5 for the Kumamoto Andisol and p = 6.17 for Hamaoka dune sand, with the value of p for the FS model being an optimized value derived from a vapor condensation experiment (Sakai et al., 2009).
For the Kumamoto Andisol (Figure 1), the agreement of the fitted WRF curves and the data was excellent for all the examined models, with the dual-(BC, VG, KO) and the VG 1 BC 2 models having R 2 values higher than .99 (Table 1). Andisols are known to be highly aggregated, with different flow processes occurring in the inter-aggregate and intra-aggregate pore systems. All of the examined models expressed this bimodality of the pore system through their first and second subfunctions. We note that the H 1 values for the Andisol were slightly greater than those for the dune sand, reflecting the size of the inter-aggregate pores. The shapes of the estimated HCF curves of those models were also very similar. Some discrepancies in the HCF curves are obvious within the transition zone (from 10 2 to 10 4 cm) because of differences in the characteristics of each function (notably between the BC function and the VG and KO functions).
Since Hamaoka dune sand has a single inflection point in the WRF (i.e., having a unimodal [single peak] pore-size distribution), the CH assumption was used for the optimized dual model in Figure 2. We confirm that ℎ b 2 or α −1 2 for the second submodel is insensitive if the parameter value lies between the air entry value and the inflection points, such as 10 or 100 cm in Figure 2. For all dual CH models, the shapes of the estimated HCF curves are similar within the range where the retention data are available (h < 10 4 cm). We note that the dual-VG-CH model is almost identical to the VG 1 BC 2 -CH model since the VG subfunction closely approximates to the BC function at large αh values. Sakai et al. (2009) showed that the FS model (designed for capillary flow using VG type retention in the wet range and for film flow using BC type retention in the dry range) described water and vapor movement with condensation and evaporation well for the Hamaoka dune sand when p was 6.17. Since the FS model also used Mualem's conductivity option of Equation 3, the VG 1 BC 2 -CH model is almost identical to the FS model as shown in Figure 2, thus validating our evaluation. Since our alternative expressions are much simpler forms of the general HCF equations, they would be very useful to describe water flow for practical applications involving a wide range of pressure heads. Furthermore, in general, the multimodal WRF with θ r = 0, asymptotically approaching θ = 0 at h is infinity, could give practically identical hydraulic conductivities in dry range as in the F I G U R E 1 Soil water retention data of Kumamoto Andisol, fitted with the dual-Brooks and Corey (BC), dual-van Genuchten (VG), dual-Kosugi (KO), and VG 1 BC 2 models (top), and relative hydraulic conductivity curves estimated according to Mualem's model with p = 0.5 (bottom) case of the WRFs with an anchoring head, ℎ m , such as the FS model.
Selection of the BC model as a second subfunction is also consistent with the Peters model where WRF was similar to the FS model assumption. Peters (2013) defined the slope of the log h − log K r curve for film flow with a parameter a, and showed that = − 1.5 based on theory of Tokunaga (2009) is reasonable for many soils. In the dual-BC-CH and the VG 1 BC 2 -CH model used in Figure 2, the corresponding slope can be approximated from the slope of the second BC function since the contribution of the first subfunction is negligible: In case of the VG 1 BC 2 -CH with p = 6.17 in Figure 2, a is −3.29. Since Mualem's definition of (q, r) = (1, 2) always results in a < −2, other sets of q and r might be necessary to have a > −2 such as (q, r) = (1, 1). In case of r = 1, the F I G U R E 2 Soil water retention data of Hamaoka dune sand, fitted with the dual-Brooks and Corey (BC)-common head (CH), dual-van Genuchten (VG)-CH, VG 1 BC 2 -CH and Fayer and Simmons (FS) models (top), and relative hydraulic conductivity curves estimated according to Mualem's model with p = 6.17 (bottom) linearly superposed conductivity function of Equation 19 can be used, similarly to the model of Peters (2013), thus allowing independent flexibility for the slope of the capillary and film conductivity curves using p 1 and p 2 . The HCF for film flow, however, needs further investigations from theoretical and experimental points of view.

CONCLUSIONS
In this study, we derived closed-form expressions based on the generalized Mualem hydraulic conductivity formulation for linearly superposed subretention functions assuming arbitrary sets of BC, VG, and KO water retention models. Because of their mathematically simple and consistent form, the derived closed-form equations are easy to implement into numerical codes to simulate water flow. We showed that different sets of sub-models can give almost identical HCFs within the range of available data. The conductivity equations can be simplified when the submodel parameters, ℎ b for the BC model, α −1 for the VG model, and ℎ m for the KO model have the same value. A hypothetical air-entry head near saturation can be introduced for multi-VG and multi-KO models to stabilize flow simulations for very fine-textured soils (having n values approaching 1). Furthermore, the superposition of conductivity subfunctions becomes linear in case the parameter r is unity.
Since the multimodal representation is very flexible to describe WRF and HCF data, the derived equations would be useful for macroporous or aggregated soils over a wide range of the pressure heads from capillary water to adsorption water. We conclude that it is important to select appropriate sets of the hydraulic models to properly analyze soil water movement.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.