A new framework to determine general multimodal soil water characteristic curves

A soil water characteristic curve (SWCC) model named as discrete-continuous multimodal van Genuchten model with a convenient parameter calibration method is developed to describe the relationship between soil suction and the water content of a soil with complex pore structure. The modality number N of the SWCC in the proposed model can be any positive integer (the so-called multimodal or N-modal SWCC). A unique set of parameters is determined by combining curve fitting and a graphical method based on the shape features of the SWCC in the log s–log Se plane. In addition, a modality number reduction method is proposed to minimize the number of parameters and simplify the form of SWCC function. The proposed model is validated using a set of bimodal and trimodal SWCC measurements from different soils, which yield a strong consistency between the fitted curves and the measured SWCC data. The uniqueness in the set of parameters provides the possibility to further improve the proposed model by correlating the parameters to soil properties and state parameters.


Introduction
The soil water characteristic curve (SWCC) describes the relationship between soil suction and water content (e.g., volumetric water content h, gravity water content w or degree of saturation S r ) of a soil. In unsaturated soil mechanics, SWCC predominates the hydro-mechanical coupling of unsaturated soils [14,25,39], since mechanical properties like the shear modulus, compression index, and yielding stress are often related to suction [1,3,31,38,45] or degree of saturation [27,50,51]. Additionally, soil properties, which are time-consuming to determine, like the unsaturated hydraulic conductivity and pore size distribution, can be derived from SWCC [2,16,26,33]. Thus, a precise description for SWCC of soils is significant for geotechnical and geo-environmental engineering, soil science as well as agriculture engineering.
A number of empirical models have been developed to reproduce the unimodal SWCC, for example, Brooks and Corey Model (BCM) [6], van Genuchten Model (VGM) [16], as well as Fredlund and Xing Model (FXM) [13]. Parameters of these models are usually obtained by bestfitting SWCC data or obtained indirectly from soil properties by using the so-called pedotransfer functions [5,19,22,36,42,43].
In recent studies, two or more pore series, resulting from the gap-graded grain size distribution or the aggregation of fine particles, have been widely observed in undisturbed soils [28], mixed soils [7,8,35] and compacted finegrained soils [10,32]. The SWCC of such soils can thus be bimodal or even multimodal, which cannot be appropriately described by unimodal SWCC models. Therefore, a set of bimodal SWCC models, as a piecewise function [8,40,47] or a continuous function [9,11,12,23,28,34,46], have been developed for the soils with heterogeneous pore structure.
The first piecewise bimodal SWCC model was developed by Smettem, Kirkby [40], who introduced two independent closed-form analytical solutions to describe the macro-and microporosity of an aggregated loam. Afterward, Wilson et al. [47] extended the method to model the hydraulic properties of a soil with three pore families by using its SWCC data. Following the framework suggested by Smettem, Kirkby [40], Burger, Shackelford [8] proposed a piecewise function to describe the bimodal SWCC of a pelletized diatomaceous earth. In their work, bimodal SWCC was divided into a macro-and a micro-sub-curve using a chosen delimiting point, leading to a piecewise bimodal SWCC function as S r ¼ S r;j þ S r;max À S r;j À Á S r;1 ; for s\s j S r;res þ S r;j À S r;res À Á S r;2 ; for s ! s j & ð1Þ where S r represents the degree of saturation; S r;max is the maximum degree of saturation; S r;res is the degree of saturation at residual state; s j is the delimiting suction; S r;j is the degree of saturation at delimiting point; S r;1 and S r;2 represent the independent closed-form analytical solutions for the local degree of saturation in the macro-and microporosity, respectively, which can be described by a unimodal SWCC function (e.g., BCM, VGM or FXM). After choosing an appropriate delimiting point, the parameters are obtained by fitting the individual unimodal function for subporosity to their corresponding SWCC data. However, in spite of the convenience in the parameter calibration process, the discontinuity feature is not expected for the numerical applications and the incorporation into a constitutive modeling. A general framework for a continuous bimodal SWCC was proposed by Ross, Smettem [34] using 'volumetric fraction approach' [21,46]. The overall pore space of the soil is regarded as the superposition of two overlapping subporosities, i.e., the micro-and macroporosity, and the bimodal SWCC function in terms of effective degree of saturation S e can be expressed as Herein, R 1 and R 2 represent the volumetric fraction of micro-and macroporosity, respectively; S r;1 and S r;2 represent the sub-SWCC curves for the macro-and microporosity, respectively, which can be described by a unimodal SWCC function. Based on this framework, a number of continuous bimodal SWCC model have been developed [9,11,12,28,34]. Details of these continuous bimodal SWCC models are summarized in Table 1. In comparison with the piecewise form, continuous bimodal SWCC model is more convenient for numerical and practical applications, but the parameters involved in these models are highly correlated, resulting in difficulties in the calibration process [18,46]. Durner [12] pointed out that the parameters should be regarded as curve shape coefficients instead of parameters with physical meanings. To best fit the parameters, a specific curve fitting procedure associated with a proper initial approximation of the parameters and appropriate constraint conditions are required. Nevertheless, identical SWCC may be reproduced by different sets of parameters [18,46], due to the possible intercorrelations among the parameters.
To overcome the difficulties in parameter calibration, another type of continuous bimodal SWCC functions has been developed based on the independent parameters related to the SWCC shape features, referred to as 'unique parameter approach' [21,46]. Gitirana and Fredlund [18] presented a bimodal model with the parameters determined with bending points from the shape of SWCC, which requires only one additional curve fitting parameter. Using a similar approach, Li et al. [23] proposed another empirical bimodal SWCC function that directly incorporates the suction and gravity water content of bending points of SWCC into the model. Wijaya and Leong [46] decomposed SWCCs into several linear segments and employed the Heaviside function to smooth the junctions of the linear segments. In this manner, all the parameters involved can be graphically determined without curve fitting procedure. Details of these models developed with the 'unique parameter approach' are summarized in Table 1. The major advantage of this type of bimodal SWCC model is the direct graphical determination of the parameters, which enables sensitivity analyses of SWCC parameters, making it possible to extend the SWCC model by relating the parameters to other soil properties and state parameters [18,46]. However, this type of model usually requires a rather complex mathematical form for the SWCC function, as shown in Table 1.
For continuous bimodal SWCC, the simple volumetric fraction models show difficulty in parameter determination, while the unique parameter approach requires a rather complex SWCC function. This problem is more significant for multimodal SWCCs. Recent studies have reported that the pore structures in unsaturated fine-grained soils are more complex than to be represented by a bimodal pore size distribution [24,30,41,44], resulting in multimodal SWCCs. To describe multimodal SWCCs, the existing bimodal SWCC models can be extended to N-modal SWCC functions. However, the parameter determination procedure is more complicated due to the significant increase in the number of parameters and the intercorrelations among them. So far, however, a general continuous N-modal SWCC model ('general' means that the modality number N can be any positive integer) is still lacking. The objective of this paper is to develop a N-modal SWCC model in a simple mathematical form with a convenient parameter calibration method.

A piecewise multimodal SWCC model-DMVGM
In this section, the piecewise bimodal SWCC function proposed by Burger, Shackelford [8] is extended to describe the multimodal SWCC. As shown in Fig. 1, the whole pore space V void of a multimodal soil is assumed to contain a permanent saturated part V res , a permanent dry part V dry , and an unsaturated part V unsat : Herein, V dry represents the volume of isolated pores in the soil, which are not accessible from outside. V res represents the volume occupied by the adsorbed water, which is strongly bonded on the soil particle surface. Thus, V res is fully saturated even at a very high suction level. V unsat represents the volume of the interconnected pores, which show varying degrees of saturation during wetting and drying cycles. From Eq. (3), the maximum degree of saturation S r,max can be determined as and the residual degree of saturation S r,res is expressed as Assuming that the unsaturated volume V unsat consists of N subporosities volumetric fraction R i for the ith subporosity can be defined as Dividing the both sides of Eq. (6) by V unsat gives in parameter calibration procedure can be found in [18] This paper Volumetric fraction approach/graphical method and curve fitting procedure N: modality number Acta Geotechnica (2021) 16:3187-3208 3191 The concept of Burger, Shackelford [8] for bimodal soils can be extended to establish the S r -s relationship for a multimodal soil by assuming an idealized N-modal pore structure as following (see Fig. 2a). For a N-modal soil, we may introduce N-1 delimiting points s i (2 B i B N) to divide the SWCC into N subcurves. One subcurve represents one subporosity. In addition, two additional 'delimiting' suctions s 1 = 0 and s N?1 = 1 are employed for mathematical convenience. For an imposed suction s 2 s i ; s iþ1 ½ Þ, it is assumed that the 1st to the (i-1)th subporosity are completely desaturated, whereas the (i ? 1)th to the Nth subporosity are fully saturated. The ith subporosity is under desaturation process, and the water volume in the ith subporosity V w,i is expressed as where S r;i represents the local degree of saturation of the ith subporosity and is described by a unimodal SWCC model (e.g., VGM). From Eqs. (3)-(9), the total water volume in the soil at a suction level s 2 s i ; s iþ1 ½ Þ, can be expressed as The degree of saturation gives Substituting Eqs. (4) and (5) into (11) gives From Eq. (12), it is not hard to find that SWCC in terms of the effective degree of saturation S e can be expressed as In this paper, VGM (with the constraint m ? 1/n = 1 suggested by Van Genuchten and Nielsen [17]) is adopted to describe S r;i : where a i and m i are the VGM parameters for the ith subporosity. From Eqs. (12)- (14), we get the degree of saturation of a N-modal soil in a piecewise form and the effective degree of saturation and the right limit we get the following inequality Therefore, DMVGM predicts a point of discontinuity at each delimiting points, as shown in Fig. 2a. For a bimodal soil (s 1 = 0, s 2 = s j , s 3 = 1), Eq. (16) degrades to and the volumetric fraction of macro (R 1 )-and micro (R 2 )porosity can be expressed as R 1 ¼ S r;max À S r;j S r;max À S r;res and R 2 ¼ S r;j À S r;res S r;max À S r;res ð21Þ Substituting Eqs. (21) in (20), it is not hard to see that the Burger, Shackelford [8] bimodal model (using VGM to describe S r,i ) is equivalent to DMVGM with N = 2, as shown in Fig. 2b.

A continuous multimodal SWCC model-CMVGM
In order to derive a continuous multimodal SWCC model, the whole pore space of the soil is regarded as a superposition of a set of overlapping subporosities, each of which occupies a volumetric fraction R i [12,34,48]. The continuous multimodal SWCC model in terms of effective degree of saturation S e is obtained by extending the bimodal function (Eq. 2) as where again N is the modality number, R i is the volumetric fraction of each subporosity with P N 1 R i ¼ 1, and S r;i is the local degree of saturation of a subporosity. This general concept was first proposed by Ross and Smettem [34]. They pointed out that S r;i may be described by any unimodal SWCC model (e.g., BCM and VGM). For the aim of simplicity, S r;i is replaced by VGM (with the constraint m ? 1/n = 1) in this study. A continuous N-modal SWCC model is then expressed as where again a i and m i are the VGM parameters for each subporosity. The SWCC function in Eq. (23) is named as continuous multimodal van Genuchten model (CMVGM).   (23)] are introduced. For practical applications of these models, convenient parameter calibration methods are required. In comparison with CMVGM, the parameter determination method for DMVGM is relatively simple due to the independence of the parameters for each subporosity, but the numerical implementation is inconvenient since DMVGM generates N-1 discontinuity points in the SWCC. Furthermore, the discontinuity feature is more significant with increasing modality number N. In contrast, CMVGM describes a continuous SWCC with convenient numerical implementation. However, including the S r,max and S r,res , the totally 3 N ? 2 parameters can be hardly determined through solely a best fitting procedure due to the strong intercorrelation among the parameters. Note that the CMVGM and DMVGM possess identical set of parameters except for the additional delimiting suctions s i for DMVGM, it is possible to calibrate the common parameters in DMVGM and use CMVGM to describe a continuous multimodal SWCC. This framework, utilizing both advantages of DMVGM (convenient parameter determination method) and CMVGM (simple and continuous mathematical function), is named as D-CMVGM. The detailed properties in D-CMVGM are shown in Table 1. In this section, the feasibility of the D-CMVGM framework for bimodal SWCCs is validated in Sect. 3.1; the determination procedure for the modality number N, delimiting suctions for DMVGM, as well as the common 3 N ? 2 parameters for DMVGM and CMVGM, is shown in Sect. 3.2; an example for reproducing a multimodal SWCC of a silty sand by D-CMVGM is demonstrated in Sect. 3.3; the development of a modality number reduction method (MNRM) is shown in Sect. 3.4.

Feasibility of D-CMVGM framework for bimodal SWCCs
The feasibility of D-CMVGM framework is validated by using CMVGM and DMVGM with identical parameters to simulate the same bimodal SWCCs, as shown in Fig. 3. Additionally, for DMVGM, the point at S e = R 2 is chosen as the delimiting point. For the silty sand [49] in Fig. 3a with a relative high ratio of a 1 /a 2 = 77 (distinct bimodal pore structure), DMVGM and CMVGM reproduce almost identical SWCC. For the undisturbed loams [28] in Fig. 3b, c with relative low ratios of a 1 /a 2 = 46 and 24, respectively, as well as the coarse sand [37] in Fig. 3d with an extreme low ratio of a 1 /a 2 = 3.6, the SWCCs reproduced by CMVGM and DMVGM are in good agreement with a slight discrepancy near the delimiting point. In general, DMVGM and CMVGM with the same set of parameters reproduce almost identical bimodal SWCC despite a slight discrepancy in a small range near delimiting point. The more pronounced the bimodal feature is, the less remarkable is the discrepancy. This phenomenon revealed a crucial fact that the parameters in CMVGM can be obtained with DMVGM based on SWCC data. Substituting the determined parameters in CMVGM then gives a continuous SWCC over the entire suction range. In fact, the parameter set of CMVGM for a multimodal SWCC may be not unique [18,46], which will be discussed in Sect. 5. The unique parameter set determined by DMVGM is regarded as one of the suitable parameter sets of CMVGM. This common parameter set of DMVGM and CMVGM can be conveniently determined in the proposed D-CMVGM framework, which is introduced in the following subsections.

Determination of S max and S res
Under the framework of D-CMVGM, the effective degree of saturation S e of a soil with complex pore structure is described by Eq. (23), while SWCC is usually represented in terms of gravity water content w, volumetric water content h or degree of saturation S r . Thus, a complete multimodal SWCC model can be generally expressed as where S is defined as general water content representing w, h or S r ; S max and S res are the maximum and residual values of that general water content, respectively. The maximum value of the water content measured during SWCC tests is adopted for the parameter S max . The parameter S res representing the residual water content at high suction range is set equal to zero. From Eq. (24), we obtain The SWCC data are then represented in terms of S e by using Eq. (25) and replotted in the log s-log S e plane to determine the other parameters. Figure 4a shows a set of unimodal and multimodal SWCCs divided into several linear segments (slope not equal to zero) and horizontal segments in the logs-logS e plane (In this paper, log X represents the base 10 logarithm of X). Under the framework of D-CMVGM, a linear segment in the logs-logS e plane is regarded as a 'subporosity', i.e., the modality number N is identical with the number of the linear segments (horizontal segments occupy zero volumetric fraction). As shown in Fig. 4a, the SWCC of silty loam [6] is unimodal; the SWCCs of kaolin-sand mixture [35], coarse sand [37] and undisturbed loam [28] are bimodal. Particularly, the SWCC of silty sand with gravel [49] is trimodal. The first linear segment represents the macroporosity, and the third linear segment represents the microporosity within the aggregations of fine particles. The second 'transition' linear segment, which is determined by the pore space in the overlapping range of macro-and microporosity, can be regarded as an extra porosity, although it occupies a small volumetric fraction. A similar finding has also been reported by Lloret and Villar [24]. They treated the microstructure of the heavily compacted 'FEBEX' bentonite as an assemblage of two distinguished porosities (macro-and microporosity) and an extra porosity in their overlapping range.

Determination of N and R i
In the logs-logS e plane, the cross-points of each adjacent linear segments and horizontal segments are chosen as delimiting points, and the volumetric fraction R i for each 'subporosity' is graphically determined. Based on this concept, a general procedure to separate a N-modal SWCC into N linear segments in the logs-logS e plane is demonstrated in Fig. 4b.

Determination of m and a for unimodal case (N = 1)
For the case of N = 1, both DMVGM and CMVGM degrade to VGM, which means that D-CMVGM with Fig. 3 Fitted bimodal SWCCs in terms of S e by using DMVGM and CMVGM with identical parameters a SWCC of a silty sand with S r,max = 0.92, S r,res = 0, data from [49], b SWCC of undisturbed Neuenkirchen loam (at a depth of 15 m) with w max = 0.46 (maximum gravity water content) and w res = 0 (residual gravity water content), data from [28], c SWCC of undisturbed Neuenkirchen loam (at a depth of 60 m) with w max = 0.42 and w res = 0, data from [28], d SWCC of a coarse sand with h max = 0.32 (maximum volumetric water content) and h res = 0 (residual volumetric water content), data from [37] N = 1 is equivalent to VGM. The SWCC in terms of effective degree of saturation S e is The slope k of the SWCC in the logs-logS e plane is defined as From Eq. (27), it is not hard to see that the slope k monotonically increases with increasing suction. Thus, the maximum slope k max is reached, when suction trends to infinite: Taking the logarithm of both sides of Eq. (26) gives When suction trends to infinite, the third term on the right-hand side of Eq. (29) vanishes, indicating a linear asymptote of VGM in the logs-logS e plane as where S Ã e is the value of effective degree of saturation on the asymptote. When the soil is saturated, the suction value s ae on the asymptote is solved from Eq. (30) by setting S Ã e ¼ 1, which gives This suction s ae is usually regarded as the air entry value of the soil. These features of VGM expressed in Eqs. (26)- (31) can also be found in [15,20].
The evolution of VGM in the logs-logS e plane is schematically demonstrated in Fig. 5a. When suction exceeds the air entry value, VGM rapidly trends to its asymptote and the slope k increases up to m/(1-m). That means, for a unimodal SWCC, the linear asymptote can be approximated by using the measured SWCC data and the parameter m is back calculated as where k is the slope of the approximated asymptote. The parameter a is the inverse of s ae (Eq. 31), which can be graphically determined, as shown in Fig. 5a. An example is shown in Fig. 5b, the SWCC of a compacted silty sand (data from [29]) is accurately reproduced by VGM with the parameters determined by the proposed calibration method.

Determination of m i and a i for multimodal case (N ‡ 2)
For the case of N C 2, the N-modal SWCC presented in the logs-logS e plane is divided into N linear segments, as shown in Fig. 4. The parameters m i and a i are determined based on the slope and position of each linear segment. In DMVGM, the slope k i of the ith subcurve in the log s-log S e plane is defined as Taking the derivative with respective to s in Eq. (16) gives Combining Eqs. (13), (14), (33) and (34), we obtain Details in the derivation of Eq. (35) are shown in Appendix A. Note that Equation (35) may be rewritten as Let us define the effective volumetric fraction R eff,i for the ith subporosity as Substituting Eqs. (38) in (37) gives For the ith subporosity, it shall be noted that the parameters m i and R eff,i are two constants. The parameter m i is characterized by the pore size distribution of the ith subporosity. That means, the slope k i solely depends on the local degree of saturation S r,i according to Eq. (39). As suction increases from s i to s i?1 , S r,i gradually decreases from 1 (fully saturated ith subporosity) to 0 (fully desaturated ith subporosity). From mathematical point of view, Eq. (39) predicts a maximum for the slope k i during the desaturation process. Furthermore, it is proved that the maximum slope k i,max exists for any combination of R eff,i (0 \ R eff,i B 1) and m i (0 \ m i \ 1) (see Appendix B). An example of the k i evolutions for different combinations of m i and R eff,i is shown in Fig. 6a (m i = 0.8) and b (R eff,i-= 0.8). Thus, the maximum slope k i,max can be generally expressed as a function of m i and R eff,i : The analytical expression for the function f m i ; R eff;i À Á is difficult to determine, but the evolution of k i,max , depending on R eff,i and m i , can be numerically obtained according to Eq. (39), which are plotted in Fig. 7a (for 0 B k The mean slope k i,mean , which is graphically determined from the ith linear segment in the logs-log S e plane, is approximated by the maximum slope k i,max for each subporosity (such an approximation is proved to be adequate to describe the SWCC accurately, see the subsequent sections): Consequently, m i can be back calculated by Eq. (41) using the prior determined mean slope k i,mean and the effective volumetric fraction R eff,i for each subporosity. Equivalently, the parameter m i can be directly determined In order to determine a i , CMVGM with prior determined R i and m i is fitted to the overall SWCC data. As shown Fig. 4b, the initial approximation for a 1 (denoted as a Ã 1 ) is regarded as the inverse of the suction of the point on the first linear segment by S e = 1. For i C 2, the inverse of delimiting suctions is adopted as the initial approximation for a i (denoted as a Ã i ). It is worth to emphasize that such an approximation is already close to a i . For this reason, the best fitting procedure is stable and converges rapidly. In this paper, the parameter a i (in the unit of kPa -1 ) is presented as the inverse of a suction value.

Summary for the parameterization of SWCCs by using the D-CMVGM framework
A general procedure to reproduce unimodal or multimodal SWCCs by using the D-CMVGM framework is summarized as following: 1. Determine S max and S res based on the measured SWCC data and calculate the effective degree of saturation S e 2. Plot the SWCC data in the logs-log S e plane 3. Divide the SWCC data into N linear segments and determine the delimiting points 4. Calculate the volumetric fraction R i and effective volumetric fraction R eff,i of each subporosity 5. Measure the mean slope k i,mean of each linear segment 6. Determine the parameters m i with prior determined k i,mean and R eff,i from the diagram in Fig. 7. 7. Determine a Ã i from the SWCC in the log s-log S e plane. Adopt a Ã i as initial approximation for a i . and fit CMVGM to all of the measured SWCC data with prior determined R i and m i .

Use CMVGM with the determined parameters to
reproduce SWCC For the case of N = 1, D-CMVGM is equivalent to VGM, and Eq. (40) degrades to Eq. (28), which corresponds to the curve for R eff,i = 1 in Fig. 7. The proposed parameter calibration method is automatically adapted to that for VGM, and hence the above described procedure is also valid for unimodal SWCCs.

Example of the reproduction of the SWCC for a silty sand with a trimodal function
The fitted SWCC of a silty sand (SW-SM with gravel according to [49]) by using D-CMVGM framework is demonstrated in Fig. 8. The maximum and residual degree of saturation S r,max and S r,res are set equal to 0.92 and 0, respectively (see Fig. 8a). Replotting the data in the logS -logS e plane, the SWCC shows a pronounced multimodal characterization. Dividing the SWCC into three linear segments (see Fig. 8b), denoted as S1, S2 and S3, respectively, a trimodal function (N = 3) is adopted to reproduce the SWCC. Setting the delimiting points at the two crosspoints of the linear segments (see Fig. 8c), the volumetric fraction R i of each subporosity is graphically determined (R 1 = 0.47; R 2 = 0.07; R 3 = 0.46) and the effective volumetric fraction R eff,i is calculated (R eff,1 = 0.47; R eff,2-= 0.13; R eff,3 = 1.0). Measuring the mean slope k i,mean of each linear segment (k 1,mean = 0.62; k 2,mean = 0.04; k 3,mean = 0.19), the parameters m i are directly obtained from Fig. 7 (m 1 = 0.78; m 2 = 0.42; m 3 = 0.17, see Fig. 8d, e). Eventually, the a Ã i (initial approximation for a i ) are directly graphically determined, and the parameters a i are obtained by fitting CMVGM to the overall SWCC data (a 1 = 1/0.12; a 2 = 1/0.50; a 3 = 1/15.5). As shown in

Modality number reduction method (MNRM) for D-CMVGM framework
In the previous subsections, the procedure of reproducing a multimodal SWCC by using D-CMVGM framework has been introduced. However, the modality number N can be further reduced by regarding the subporosity with relatively low volumetric fraction as the overlap of its adjacent subporosities. Such a subporosity is named as 'transition subporosity', and its linear segment is named as 'transition segment (TS)'. In this study, the subporosity, whose volumetric fraction is lower than 0.1, is regarded as the transition subporosity. The middle point of a transition segment is chosen as a delimiting point to split the transition segment. By incorporating half of the volume into the former adjacent subporosity and the other half into the latter adjacent subporosity, the modality number is reduced (modality number coincides with the number of linear segments). This method is named as modality number reduction method (MNRM), which simplifies the form of SWCC function and reduces the number of unknown parameters. A schematic representation for dividing a N-modal SWCC into linear segments and transition segments in the log s-log S e plane by using MNRM is shown in Fig. 9.
The procedure for reproducing the SWCC of a silty sand (SW-SM with gravel according to [49]) by using D-CMVGM framework incorporating MNRM is shown in Fig. 10. The second subporosity with a volumetric fraction of 0.07 is regarded as a transition subporosity and its linear segment as transition segment. The corresponding parameter calibration procedure is shown in Fig. 10. As can be seen, the SWCC is precisely reproduced, while the SWCC function is simpler, and less parameters are required.

Simulation of bimodal and trimodal SWCCs of mixed soils
In Fig. 11, the measured SWCCs of four mixed soils (S1-S4) from [35] are demonstrated. The four artificial soils were prepared by mixing coarse kaolin and Ottawa sand with different fines contents, and therefore they have different pore structures and different SWCCs. As shown in Fig. 11a, the SWCC of S1 can be regarded as an assembly of three linear segments. The second subporosity occupies a volumetric fraction of 0.08, which is less than 0.1. However, in order to validate the ability of the D-CMVGM to reproduce a trimodal SWCC and improve the accuracy of the curve fitting, this subporosity is not seen as a transition subporosity, and a trimodal function is used to reproduce the SWCC of S1. In Fig. 11b, the SWCC of S2 splits into four linear segments. As the second linear segment is almost horizontal (the volumetric fraction of second subporosity is close to zero), MNRM is applied, and the middle point of the second linear segment is chosen as a delimiting point. Therefore, a trimodal function is used to describe the SWCC of S2. The SWCCs of S3 and S4 splits into three linear segments, as shown in Fig. 11c, d, respectively. For the same reason as Fig. 11b, the second linear segments of both SWCCs are regarded as a transition segment, and delimiting points are set in the middle of the transition segments, i.e., a bimodal function is adequate to reproduce the SWCCs of S3 and S4. The volumetric fraction R i , the mean slope k i,mean as well as a Ã i (as the initial approximation of a i ) for each subporosity are presented in Fig. 11. The parameters for D-CMVGM framework as well as the h max and h res for each sample are listed in Table 2. Figure 12a, b demonstrate a good consistency between the fitted curves and the measured SWCCs.

Simulation of bimodal SWCCs of a silty sand
Angerer [4] prepared a set of reconstituted samples of low plasticity silty sand, which were statically compacted to different initial densities (I d = 0.5, 0.7 and 0.9) at different water contents (w = 3%, 6% and 10%). The fine content of the soil is 9.5%, including 1% clay and 8.5% silt. The SWCCs of the samples were measured over a wide suction range up to about 1 9 10 6 kPa by using both suction tensiometers (for suction lower than 1 9 10 3 kPa) and a chilled mirror hygrometer (for suction higher than 1 9 10 3 kPa). In this paper, only the SWCCs of the samples compacted at a medium density of I d = 0.7 are presented and reproduced.
In Fig. 13a, the SWCC for the soil compacted at w = 3% is demonstrated in the log s-log S e plane, which consists of two linear segments and a transition segment. A delimiting point is set in the middle of the transition segment according to MNRM. Therefore, a bimodal function is used to describe the SWCC. A similar approach is applied to the SWCCs for the samples compacted at w = 6% and w = 10% in Fig. 13b, c, respectively. The parameters of the three SWCCs are shown in Table 3. With the S r,max and S r,res , the fitted SWCCs in terms of degree of saturation along with the SWCC data are presented in a conventional log s-S r plane in Fig. 13d, by which a good agreement between the fitted curves and measured SWCCs is shown.
In addition, the influence of the compaction water content on the SWCCs can be investigated from the variations of the parameters in Table 3. It is noted that the parameters R 1 , R 2 , a 1 and a 2 are affected, while the other parameters remain almost constant. When regression analysis is conducted and parameters R 1 , R 2 , a 1 and a 2 are related to compaction water content, the SWCC for the soil compacted at other water content can be estimated.

Discussion of the uniqueness of the set of parameters
As mentioned in Sect. 3, the parameters involved in CMVGM are highly correlated, and thus the use of the least square fitting approach for parameter determination might cause convergence problems in the optimization process [18,46,52]. Gitirana Jr and Fredlund [18] pointed out that a unique set of parameters may not exist, when the fitting parameters is not related to the shape features of curves. In this work, this issue is analyzed by reproducing the identical SWCC using CMVGM with two different set of parameters.
In Fig. 14a, the SWCC of a silty sand (SM with gravel according to [23]), presented in the log s-log S e plane (w max = 0.176 and w res = 0), is approximated by three linear segments. Two delimiting points are set at the crosspoints of the linear segments. The parameters determined for a trimodal SWCC function are referred to as 'solution 1'. In order to find another set of parameters, a specific SWCC separation approach is introduced in Fig. 14b, where the second linear segment is regarded as a transition segment. One delimiting point is set in the middle of the transition segment, and the other delimiting point is used to divide the former first linear segment into two parts. Based on this specific separation approach, the determined parameters are referred to as 'solution 2'. The parameters in both solutions are shown in Table 4.
The fitted SWCCs with both sets of parameters are demonstrated in Fig. 14c (in terms of effective degree of saturation S e ) and d (in terms of gravity water content w). Despite the different parameters, the two fitted curves are almost identical and consistent with the SWCC data, revealing a crucial fact that the set of parameters for CMVGM may be not unique for identical SWCC. Potentially causing the convergence problem and uncertainties in the curve fitting procedure, this shortcoming of CMVGM is overcome under the proposed D-CMVGM framework. A unique set of parameters can be determined with a predefined SWCC linearization and separation procedure in the log s-log S e plane. This feature provides the possibility to

Conclusion
A continuous N-modal SWCC model D-CMVGM with a convenient parameter calibration method is proposed, by which the modality number N can be any positive integer. The CMVGM provides a continuous function to describe the multimodal SWCC of the soils with heterogeneous pore structure. However, the determination of all the parameters solely with a curve fitting procedure leads to convergence problems and enhanced uncertainties, due to the none-uniqueness in the parameters of CMVGM. This problem is overcome under the developed D-CMVGM framework. A unique set of parameters are conveniently determined by a prior SWCC linearization and separation procedure in the log s-log S e plane. The modality number N corresponds to the number of linear segments of the SWCC presented in the log s-log S e plane. In addition, MNRM is proposed to reduce the number of parameters and simplify the form of SWCC function. The parameters R i and m i can be graphically determined, and the parameters a i are determined using a curve-fitting procedure with known parameters R i and m i . Eventually, the parameters are substituted into CMVGM to reproduce a continuous multimodal SWCC.
The mathematical form of D-CMVGM is relatively simple in comparison with other multimodal (bimodal) Fig. 10 A fitted multimodal SWCC of a silty sand from the D-CMVGM framework incorporated with MNRM (data from [49]) SWCC model developed with unique parameter approach. In this work, a total of 9 bimodal SWCCs and 3 trimodal SWCCs from different soils are reproduced, and the fitted curves show good consistency with the SWCC data. In addition, the D-CMVGM framework can be further Table 2 Parameters of D-CMVGM for soils S1-S4 Applying VGM for the ith subporosity (Eq. 14) to Eq. (42), we obtain Appendix B: proof for the existence of the maximum slope k i,max a) For the case 0 \ R eff,i \ 1 (i \ N) In Eq. (39), k i is continuous on S r,i [ [0,1] and differentiable on S r,i [ (0,1). After the Lagrange's mean value theorem, there is a value n of S r,i [ (0,1) such that dk i dS r;i S r;i ¼ n Thus, the slope k i at S r;i ¼ n is the maximum slope k i,max . b) For the case R eff,i = 1 (i = N or unimodal SWCC) If R eff,i = 1, Eq. (39) degrades to Fig. 12 Fitted SWCCs in terms of volumetric water content from the D-CMVGM framework (data from [35]) a for soils S1 and S2, b for soils S3 and S4 Fig. 13 a Separation of SWCC w = 3% into linear segments, b separation of SWCC w = 6% into linear segments, c separation of SWCC w = 10% into linear segments, d reproduction of the SWCCs of a medium dense silty sand compacted at different water contents (data from [4]) Thus, the slope at S r;i ¼ 0 is the maximum slope k i,max . Combining a) and b), the existence of the maximum slope k i,max for any combination of m i [ (0,1) and R eff,i-[ (0,1] is proved. Acknowledgements The support of the China Scholarship Council (Number 201608080128) is greatly acknowledged by the first Author. Availability of data and material The data used to support the findings of this study are available from the corresponding author upon request.
Code availability The code used to support the findings of this study are available from the corresponding author upon request.  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.