Calibration of BBM Parameters using the Modified State Surface Approach

. The Barcelona Basic Model (BBM) developed by Alonso et al. [1] is the first and the most widely used elasto-plastic model for unsaturated soils. The BBM successfully explained many key features of unsaturated soils and received extensive acceptance. However, there is lack of a well-established method for selecting parameter values for the BBM from laboratory tests, although a variety of methods have been recently developed for calibrating model parameters for the BBM. Concerns still exist on the correctness and robustness of such parameter value selection procedures. The above statements were evidenced by a recent benchmark exercise on selection of parameter values for the BBM organized within a "Marie Curie" Research Training Network on "Mechanics of Unsaturated Soils for Engineering" (MUSE)[2]. Experienced constitutive modelers in unsaturated soils from 7 prestigious teams were provided with the same experimental results on an unsaturated soil to calibrate the parameter values in the BBM. Theoretically, the calibrated parameters from different teams are expected to be the same or at least very close. However, the selected parameter values by the 7 teams are surprisingly widely different. This paper first discussed the limitations in the existing methods to calibrate the parameter values in the BBM. A novel approach was then proposed to calibrate the parameter values for the BBM. The approach takes advantage of the close-form solution of the BBM, which is derived based upon a newly proposed Modified State Surface Approach to study the unsaturated soils [3-6]. The same experimental data, used in the MUSE benchmark exercise, were reanalysed using the proposed approach to calibrate parameters for the BBM. The results were compared with those in the MUSE benchmark exercise from which the simplicity, effectiveness, and robustness of the proposed method were evaluated.


Introduction
The Barcelona Basic Model (BBM), first proposed by Alonso et al. [1] in 1990, has received extensive acceptance and been successfully applied to explain various key features of unsaturated soils' behaviour. It has also been implemented into Finite Element software to analyse boundary value problems related to earthworks [7], railways and climate [8,9], field tests [10], and disposal of nuclear waste underground [11]. However, dissemination and use of the BBM outside the unsaturated soils research context have been very limited. Possible contributory factors in this included, but not limited to, (1) the experimental difficulties associated with performing suction-controlled tests that are needed for calibrating the model parameters, and (2) uncertainty in how best to select BBM model parameter values from laboratory test data and concerns on the robustness of such parameter value selection procedure [12]. While some advances have been introduced to the former topic to characterize unsaturated soil behaviour in a faster and more accurate way [13][14][15], little progress was made in the latter one. Despite the fact that a great number of efforts have been recently dedicated to calibrating model parameters for the BBM [2,[12][13][14][15][16][17], concerns still exist on the correctness and robustness of such calibration procedures. The above statements were evidenced by a recent benchmark exercise on selection of parameter values for the BBM organized within a "Marie Curie" Research Training Network on "Mechanics of Unsaturated Soils for Engineering" (MUSE). 7 teams of experienced constitutive modelers provided with the same experimental results on an unsaturated soil to calibrate the parameter value for the BBM. Theoretically, the calibrated parameters were expected to be the same or very close, yet the selected values proved to be surprisingly scattered [2].
In this paper, a critical review was first performed for the existing methods to calibrate the BBM model parameter values in which the limitations of these methods were presented. It is found that all the current methods were built upon the incremental forms of elasto-plastic theory in which the original BBM was proposed. The limitations of incremental forms in providing an overall best fit to the BBM were explained and discussed.
Based on the discussions, a novel approach was proposed to calibrate the model parameters for any constitutive model including the BBM. The proposed method is based upon the Modified State Surface Approach (MSSA), which is proposed by Zhang and Lytton [3][4][5][6] to investigate elasto-plastic behaviour for unsaturated soils. Using the MSSA, the closed-form surface equations were derived for the BBM, which can help to establish the objective function for the BBM model calibration. Later, the Broydon-Fletcher-Goldfarb-Shanno (BFGS) search algorithm, which is a quasi-Newton method, used to minimize the objective function by searching for the best set of model parameters [23]. The derived closed-form objective functions using the MSSA along with the BFGS method, was then implemented to calibrate the parameter values in the BBM. The same test results used in the MUSE benchmark exercise were used to calibrate the model parameters for the BBM [2,24]. The results were compared with those in the MUSE benchmark exercise from which the simplicity, effectiveness, and robustness of the proposed method were evaluated.

BBM: An Overview
The BBM was developed by Alonso et al. [1] using mean net stress p, deviator stress q and matric suction s as stress state variables, where p is the excess of mean total stress over pore air pressure and s is the difference between pore air pressure and pore water pressure. It is an extension of the Modified Cam Clay model for saturated soils [25] into unsaturated soils. In the elastic zone, the following equation was proposed for the variation of specific volume due to small changes of the net normal tress p and suction (Please see the list of symbols for definitions).
( ) e s at dv dp ds p s p   =+ + (1) The elastic strain induced by changes in q is computed using Eq. (2): 1 3 e s d G dq   =   (2) A loading collapse (LC) and suction increase (SI) yield curves were proposed to separate the elastic and plastic zones as follows under isotropic conditions:      (4) where, (5) It is worth noting that the SI yield curve was mostly abandoned in the later refinement of the BBM [4]. As a result, any features related to the SI yield curve is not discussed in this paper for simplicity purposes.
The BBM was then extended into triaxial stress states by assuming the yield curve at constant suction s is an ellipse in the p-q plane, similar to the Modified Cam Clay model [2], as follows: where p s =ks and represents the soil cohesion increase with suction, and it was assumed to be a linear relationship. The yield surface can also be expressed as the following function by combining Eqs. (3) and (6) The BBM assumed purely volumetric hardening. The plastic deformation of the specific volume is solely dependent on the degree of expansion of the yield curve and it is independent of the actual stress path employed to produce the expansion (Wheeler and Karube 1996) (8) Under continuous shearing, the BBM assumed that the soil will ultimately attain a critical state and the soil will fail if the following condition is met: ( ) q M p ks =+ (9) A non-associated flow rule was used to calculate the deviatoric plastic strains as follows: In the original BBM, it was suggested that  is chosen in such a way that the flow rule predicts zero lateral strain for stress states corresponding to Jaky's K0 values [26] Zhang et al. [15] argued that K 0 for unsaturated soils cannot be a constant and derived the close-from solution for the K0 for unsaturated soils. Consequently Eq. (11) should not be used. Instead, it is proposed that  be used as an additional constant for the "modified" BBM.
The BBM used the following equations to calculate the total volumetric and shear strain increments: (13) Where  stands for strain, the superscripts "e" and "p" represent elastic and plastic strains respectively, while the subscripts "v" and "q" stand for volumetric and deviator strains, respectively. Under isotropic conditions, the BBM was also presented an integrated form. For isotropic compression tests at constant suction, s, the BBM adopted the following equations to predict soil behaviour at virgin states:

Existing Methods for BBM Parameter Calibration and Limitations
In constitutive modelling, the process of calibrating model parameters is based on experimental tests that encompass the range of stress states for the specific boundary value problem being examined [18][19][20]. The ultimate objective of the calibration is to determine the most appropriate parameters that can accurately predict the soil's response to the available experimental data. It is recommended to have redundant data points during the calibration process, especially when dealing with highly nonlinear soil behaviour, in order to obtain the most accurate representation of the soil's behaviour. As a result, calibration of a model inevitably requires statistical analysis. The calibration of an elasto-plastic model is a mathematical process that involves minimizing the error between experimental results and numerical predictions through the adjustment of model parameters. Physical constraints may apply to the model parameters in certain situations, thereby making the calibration process a constrained optimization problem.
The optimization procedure consists of two main steps, Firstly, the formulation of an objective function denoted as F(X), which measures the difference between the theoretical and experimental results. It is essential to include the stress/strain behaviour at every point in each test to obtain an accurate representation of the soil's overall behaviour. Secondly, the selection of an optimization strategy enabling the search for the minimum of the objective function.
(3). Use the yield loci from isotropic normal compression tests at different suction levels to fit Eq.
(4). Use several shear strengths tests at different suction levels to fit Eq. (9) to obtain M and k. The method is straightforward. Several objection functions were clearly defined, and standard optimization techniques can be used to obtain the best fit. However, as pointed out by Gallipoli et al. [12], in the BBM, different aspects of soil behaviour are each affected by more than one of these parameters while, at the same time, a single parameter controls more than one aspect of soil behaviour. This makes it difficult to associate each parameter with a single aspect of material behaviour. Because of this, iterative approaches have often been needed in which some parameters are initially fixed by matching one particular aspect of soil behaviour but are subsequently adjusted to improve prediction of other experimental features. Such iterative approaches involve a substantial degree of personal judgement, often leading to the selection of widely different parameter values from the same experimental data.
Gallipolli et al. [12] proposed an improved sequential method especially designed for calibration the 5 model parameters in the BBM, β, λ(0), r, p c , and N(0) under isotropic conditions. Degrees of freedom in the BBM are progressively eliminated in a specific order, so that the corresponding parameter values are selected one at a time without having to make assumptions about the values of remaining parameters. The proposed calibration method is simpler and less subjective than conventional procedures. Due to the limited need for personal judgment, the proposed procedure makes it possible for a relatively inexperienced user to select parameter values from a set of suction-controlled laboratory tests in a relatively robust and objective fashion.
Instead of using the sequential method, some researchers have opted for global optimization methods to calibrate model parameters for constitutive models that involve numerous parameters [18,19], just like the BBM. These methods involve minimizing an objective function that measures the difference between model predictions and laboratory data and are theoretically more rigorous than the sequential approach. However, if the methods are to be used for BBM model calibration, they require an in-depth understanding of every detail in the BBM to establish the objective function from multiple individual constitutive relations, such as Eqs. (1) through (13) in Section 2. Additionally, data preprocessing requires substantial personal judgment, and the use of advanced numerical methods in engineering optimizations is necessary to solve the complex objective function. These methods often operate as "black boxes," hiding the physical meaning of individual parameters and requiring significant computational effort. Even experienced users may find it challenging to use these programs, as they need to be customized specifically for calibrating model parameters in the BBM [13].
These challenges have been recognized by many researchers. A benchmarking exercise on calibration the model parameters for the BBM was organized within a "Marie Curie" Research Training Network on "Mechanics of Unsaturated Soils for Engineering which was funded by the European Commission from 2004 until 2008 [27]. Researchers from 7 prestigious universities took part in the benchmarking exercise: the University of Glasgow (GU), UK; the University of Durham (DU), UK; the Università degli Studi di Trento(UNITN), Italy, the École Nationale des Ponts et Chaussées(ENPC), France; the Università degli Studi di Napoli Federico II(UNINA), Italy; the Universität Innsbruck(UNINN), Austria; and the University of Strathclyde(USTRAT), UK. The first 5 of these were members of the MUSE Network and the last 2 were external participants. The exercise was coordinated from the University of Glasgow (GU). D'Onza et al. [2] elaborated these efforts in great detail. In practice, 6 of the 7 teams except UNINN used a sequential approach to isolate specific features of behaviour to determine the values of different individual model constants, but then employed some degree of iteration or compromise. Different teams also chose to place greater or lesser emphasis on aspects of behaviour or on individual tests. Different from the other universities, UNINN performed a formal global optimization process using inverse analysis. This involved simultaneous optimization of the values of most of the 10 soil constants under isotropic conditions and the initial value of the hardening parameter, by attempting to minimize suitable objective functions describing the differences between model simulations and experimental results. Exceptions were G, M and k, which were determined by the UNINN team in a more conventional fashion. Theoretically, the calibrated parameters from different teams are expected to be the same or at least very close. However, the selected parameter values by the 7 teams are surprisingly widely different, which is presented in Table 1. For example, p c varies from 10 -4 to 1017, with a difference of 10 7 .

A simple example: plane fitting
Before discussing the limitation of existing methods for the BBM calibration and possible solutions, let us have an exercise on fitting a simple plane to demonstrate how a model calibration is performed with different approaches. Assume that we have a set of four points of A, B, C, and D in the threedimensional (3D) space as shown in Fig. 1 and wants to best-fit them with a plane. As shown in Fig. 1, points A and B are in the xOz plane with x-, y-, and zcoordinates of (1,0,2) and (2,0,1) respectively, while points C and D are in the yOz plane with coordinated of (0, 0.5,1.5) and (0,1,0.5) respectively. Alternatively, the four points can be viewed as a group of test results and the plane can be viewed as a model for the test results with an expression as follows: The model can also be written in terms of an incremental form as follows: dz adx bdy =+ (17) Where z a x  =  (18) z b y  =  (19) Four approaches can potentially be used to calibrate the model parameters a, b, and c in Eq. (15) as follows: (1) a sequential approach based upon the incremental form as represented by Eqs. (18) and (19), (2) a global optimization approach based upon the incremental form as represented by Eq. (17), and (3) a global optimization approach based upon the closed form solution as represented by Eq. (16).
The details of the three different calibration approaches are described as follows.
For approach (1), one can use the results of points A and B in the xOz plane to find the value of the model parameter a. The objective function is: Since points A and B lie on a line of z=3.00-1.00x and there are only two points, the slope of the line is / 1.00 z x a   = = − (21) Similarly, one can use the results of points C and D in the yOz plane to find the value of the model parameter b. Since points C and D lie on a line of z=2.500-2.00y, the slope of the line is  Table 2 shows a summary of the results obtained from three different approaches. The R-squares (coefficient of determination) for the two fitting processes are 1.00 since only two points are used to fit a line in both cases. c obtained from the two regressions were 3.00 and 2.50, respectively. One inconsistency can be noticed here that two different values for the parameter c are calibrated. When the two c values are used, the predicted z values both have overall R-squares of 0.60, which are significantly lower than the 1.00 of R-squares obtained from the two individual optimizations. For approach 2, the objective function to best fit Eq. (16) by using the four points is as follows: Where (dz) mi and (dz) pi is the measured and predicted results of dz for a point relative to the initial condition, respectively. In order to perform the global optimization, one test point must be selected as the initial values for calculating dx, dy, and dz to perform the least square optimization to obtain the parameters a and b. It is worth noting that for different initial conditions, the results obtained will be different. For example, as shown in Table 2, if point A is selected as an initial point (option 1), a and b obtained from the optimization process using Eq. (23) are -0.9167 and -2.5000, respectively. The R-square for the fitting is 0.9167. The corresponding c is 2.9167 based upon point A. The overall coefficient of correlation is 0.9667. If points B, or C, or D is used as the initial conditions, the R-squares for the global optimization change from 0.9444 to 0.9762, and the overall Rsquares vary from 0.9667 to 0.9778. In other words, the results depend upon the selection of the initial conditions. Of course, the easiest and most accurate way is to use the coordinates of the four points and the closed-form as represented in Eq. (15) to perform a standard leastsquare analysis (approach 3). The object function for the fitting assignment can be derived as follows: Where z mi and z pi is the measured and predicted results of z based upon the model parameters a, b, and c at the different values of x and y, respectively. Without providing more details, the results of the optimization are a=-0.8500, b=-2.3000, and c= 2.7500 with a plane has the following expression: 2.7500 0.8500 2.3000 The overall coefficient of correlation for the fitting is 0.9800, which is the highest value in all three approaches as shown in Table 1. Obviously, approach (3) is the correct and simplest approach for the calibration process.
The above example is simple, but it clearly demonstrates the limitations in the first two approaches. Approach (1), a sequential method, used partial results to perform calibrations for part of the model. Two independent optimization processes are performed to get the model parameters, and results are then combined for the whole model. It can be seen from the above example that each optimization has high R-square and best-fit locally, but the results do not represent the overall best-fit for the whole model. There is an inconsistency with the values of c as well. Approach (2) is a global optimization process to best fit the whole model based upon the incremental formulation. Application of incremental formulation requires specifying an initial point, and the results are highly dependent on the selection of the initial points. This is attributed to the fact that this approach implicitly assumes that the initial condition is perfectly on the plane. However, as can be seen from Eq. (24) and Fig. 1, none of the four points is exactly fallen on the plane. The calibration results using approach (2) can assure 100% fit of the initial point and a best fit of the derivatives of Eq. (16) relative to the initial point. However, unless the initial point is exactly fallen on the plane, approach (2) cannot obtain the best fit of the model represented by Eq. (1). Approach (3) is the most rigorous way and can get the overall best-fit for all results, although the R-square is lower than those in some of the partial fittings.
In order to make the approach (2) work, a "modified" global optimization approach (4) based upon the incremental form as represented by Eq. (17) can be used to obtain the correct results. Instead of using a test point as initial condition, an additional parameter of predicted z value z p0 at the initial conditions of x and y must be introduced into the calibration process. The corresponding objective function can be written as follows: For example, when point A is used as an initial condition, only the independent variables x and y values of (1,0) of pint A are used as initial conditions, while the dependent variable z is considered as an unknown value for optimization. The results obtained are the same as those in Approach (3). It is worth noting that calibration based upon Eq. (25) can obtain correct a and b values no matter if the initial point is on the plane. The downside is that when the results are used to make predictions, one must start from the same initial x and y values used in the calibration to assure the equivalent c value is correct and obtain the correct prediction results. This is easy to do for this simple plan-fitting but will become quite difficult to achieve when dealing with a much more complicated model like the BBM with many tests (initial conditions) as well as elasto-plastic behaviour involved. More discussions can be provided in later sections in this regard.

Limitations in the MUSE Benchmark Exercise
From the above discussions, when examining the MUSE benchmark exercise [2], it is found that the following limitations exist in the existing methods: 1. All the researchers used the sequential approach including the UNINN team in the MUSE benchmark exercise. Although each of the optimization represents best fits for specific features, the combined results do not represent an overall best fit of the whole model.
2. All the researchers used the incremental formulations when calibrating the BBM. To calculate the increments, all calculations must start from certain initial conditions which are implicitly assumed to be perfect with no errors, which is problematic. Since many tests are used to calibrate the BBM parameters, many of such initial points are used.
3. In all the calibrations, it is assumed that all the tested soil specimens have the same stress histories. This is evidenced by the fact that only one preconsolidation stress was reported by each team for all the tests. As can be seen in the later sections, the tested soil specimens most likely have different stress histories. Existing methods also require the preconsolidation stresses are determined accurately, which is difficult to achieve using the Casagrande method [36] (Please see Zhang et al. [28] for more discussion).

Modified State Surface Approach
The exercise on the plane-fitting clearly demonstrated the benefits of model calibration based upon the analytical solution of the model (plane). It also demonstrated that if the incremental formulation is used, one cannot easily perform a correct model calibration, no matter if a sequential method or global optimization method is used and what point is used as the initial condition. Because of these shortcomings, the use of simple analytical methods should remain the preferred option for selecting the values of a constitutive model.
Zhang and Lytton [3][4] developed the MSSA to explain the elastoplastic behaviour of unsaturated soils under isotropic stress conditions. The closed-form surface equations for the BBM was derived based upon the principle of the MSSA. Afterward, the close-form solutions for the BBM under triaxial loading conditions was also derived [6]. Zhang & Lytton [5] extended the MSSA for the coupled hydro-mechanical behaviour of unsaturated soils. Riad & Zhang [17,29] further extended the MSSA to include the coupled hydromechanical hysteresis for unsaturated soils. In this paper, only the relationships for volume changes in soil structure are discussed. In the following section, the principles of the MSSA are introduced, and then the method is utilized for derive the close form solutions for the BBM for the model calibration purposes.

Principle of the MSSA
The principle of the MSSA can be illustrated by Fig. 2 [3]. Fig. 2 shows the stress paths for three isotropic loading-unloading-reloading tests. At an arbitrary constant suction s = s 2 , the soil specimen initially resides at point D with an initial yield curve of LY 1 . The preconsolidation stress is p 0 * at s = 0 kPa, and the yield stress at s = s 2 is represented by point E in Fig. 2b. The soil is then loaded from D to E and subsequently to V, unloaded from V to D', and finally reloaded to F in Fig. 2b. This figure presents a typical soil response in the v-lnp plane, neglecting hysteresis. The process yields the following observations: 1. The virgin loading curve EVF's shape and position remain constant in the v-lnp plane, regardless of stress path or history. Plastic loading merely expands the curve's range, such as the initial EVF curve becoming VF after loading from D to E to V.
2. In the v-lnp plane, the unloading-reloading curve retains its shape and position during elastic loading or unloading (e.g., D to E, V to D', or D' to V). For plastic loading, the curve's shape and slope remain constant, but its position shifts downward in parallel with the original curve. The elastic zone's range expands due to the increase in preconsolidation stress from p 0 * to p 1 *.
3. The yield point is the intersection of the unloading-reloading curve and the virgin loading curve for s=s 2 , e.g. points, E, V, and F.
Figs. 2a and 2b illustrate compression tests conducted at two additional arbitrary suction levels of s 1 and s 3 , such as stress paths ABC and GHI, along with their results. Similarly, in the elasto-plastic zone, the shapes and positions of the virgin curves at different suction levels such as BC, EF, and HI as shown in Fig. 2b will always be the same if they are plotted in the v-p-s space as shown in Fig. 2c. Continuously conducting tests at varying suction levels results in a fixed "plastic (virgin) loading surface" in the v-p-s space, such as BEHUXYZWB in Fig. 2c, which is unique. The uniqueness of the state boundary surface is a fundamental assumption in the constitutive modelling of elastoplastic soil behaviour, and experimental evidence has validated the uniqueness of the state boundary surface for unsaturated soils [30]. The plastic surface BEHUXYZWB in Fig. 2c is the shape of the state boundary surface for isotropic conditions. In the meantime, as long as the stress path is in the elastic zone, the soil behaviour is stress path independent and soil responses can be represented by a unique surface. In the v-p-s space, the following assertions can be made for the elastic and plastic surfaces: 1) The shape and position of the plastic surface are always the same for the soil in the v-p-s space. Virgin loading only changes the range of the plastic surface.
2) During an elastic loading or unloading process, the shape and position of the unloading-reloading elastic surface and the plastic surface remain unchanged in the v-p-s space. The specific volume of any isotropic elastic loading or unloading stress path must fall on the elastic surface in the v-p-s space. During a plastic loading process, the shape of the unloading-reloading elastic surface remains unchanged ( and  s are constants for an assumed planar elastic surface), but its position will shift. Specifically, the unloading-reloading elastic surface will move downward in parallel with the original unloadingreloading elastic surface. For example, the unloading stress path V to D' will fall on the new elastic surface A'WVUG'D'. The volume change of any isotropic plastic loading stress path must fall on the plastic surface in the v-p-s space. 3) The yield curve is the intersection of the unloading-reloading elastic surface and the plastic surface.
The MSSA was extended to the triaxial stress states in the v-p-q-s space (which also considers deviator stress, q) as follows [6,15]: 1) There is a unique state boundary surface in the elastoplastic region which is always unchanged in the v-p-q-s space.
2) The elastic surface is movable, but only moves when there is plastic loading. The elastic surface is fixed when there is elastic loading or unloading.
3) All the soil responses will fall on either the elastic or plastic surface.
4) The intersection between the elastic and plastic hypersurface is the yield surface; and 5) The plastic hypersurface ends when the soil fails, which is at the critical state.

Analytical solution of the BBM
Using the MSSA, Zhang and Lytton [3] derived the closed-form expressions for the BBM under isotropic stress conditions. Fig. 3 shows the elastic and plastic surfaces used in the BBM. They include an elastic surface AFIH and the plastic surface which is made up of two parts: a plastic collapsible surface FIJG (corresponding to the LC yield curve) and a plastic expansive surface HIJC (corresponding to the SI curve). The elastic surface AFIH as shown in Fig. 2c has the following expression: The plastic collapsible surface FIJG can be expressed as: (25) The plastic expansive surface HIJC has the following expression: = atmospheric pressure, = slope of the unloadingreloading line associated to the mean net stress,  s = slope of the unloading-reloading line associated with soil suction, and C 1 and C 3 = are constants. The superscripts "e" represents the elastic change in the specific volume.
Zhang [6] also derived the close-form solution for the BBM under a triaxial stress states for the plastic collapsible surface as follows:  (27) where q = 1 - 3 deviatoric stress, k = parameter that relates cohesion and suction, and M = slope of theoretical critical state line. Under triaxial stress states, the elastic and plastic expansive surfaces ae the same as Eqs. (24) and (26) respectively.   (26) for the BBM model within the framework of the MSSA respectively. Partial derivatives of the surfaces with respect net normal stress p and suction are also shown in Fig. 3. These partial derivatives are consistent with the original BBM [1] as discussed in section 2. In other words, Eqs. (24) through (26) provide integrated closeform solutions for the BBM under isotropic conditions, while Alonso et al. [1] mostly used an incremental formulation to develop the BBM in which partial derivatives were first proposed to represent specific features of unsaturated soil behaviour along some special stress path, and then the partial derivatives were assembled to the whole model of BBM according to critical state soil mechanics theory to represent soil behaviour under all possible stress paths. It is worth noting that use of incremental formulation to develop a constitutive model is a well-established approach for researchers from the birth of elasto-plastic theory. However, it results in great challenges for model parameter calibration as discussed in section 3 for plane fitting, yield curve determination as discussed in Zhang et al. [28], and laboratory testing as shown in Zhang [31].

Relationship between the close-form solution and the incremental formulation
Zhang and Lytton [3] used the MSSA to successfully replicate all the examples in the Alonso et al. [1], indicating that the MSSA is consistent with the existing theories of elasto-plasticity for unsaturated soils. In addition, the MSSA clearly explains the relationship between different components of a critical state constitutive model. Especially, the advantage of having the uniqueness of the state boundary surface can be fully taken to significantly simplify the calibration process. This, combined with the closeform solutions, makes the calibration of the BBM model parameters much easier and more straightforward.

Use of the MSSA to Calibrate the BBM Parameters: Formulation of Objective Function
The exercise involving plane-fitting highlighted the challenges of calibrating models based on incremental formulation. The task becomes even more arduous for complex models like the BBM, which comprises over a dozen constitutive relations for different features of unsaturated soils. Obtaining accurate model parameters using incremental formulation becomes an overwhelming challenge in such cases.
In contrast, calibration of models through an integrated close-form solution simplifies the calibration process significantly. Close-form solutions for the BBM were obtained through the MSSA approach [3,6]. Using these solutions for BBM model calibration is much more straightforward and less confusing. The following sections discuss how to establish objective functions for BBM model calibration in different scenarios. In the following calibration process, it is assumed that the elastic and elasto-plastic behaviours of the unsaturated soil is independent of each other. Consequently, two optimization processes were performed for model parameter calibrations, one is for the elastic zone, and the other is for the elasto-platsic zone. In case that the elastic and elasto-plastic behaviours are not independent of each other, then only one optimization is needed and all test results should be included in one single calibration process.
In general, calibration of the BBM parameters is a process to find a suitable combination of model parameters of ( ) ( ) (28) Which can minimize the differences between the measured results from all the tests and the predicted results by the BBM using the model parameters as selected for Eq. (28) as follows: Where n represents number of tests performed on the unsaturated soils, w ij = weighted factor, v mi = measured specific volume, v pi = predicted specific volume, ( ) It is worth noting the preconsolidation pressure is implicitly included in Eq. (28). According to the MSSA, yield curve/surface is the intercept of the elastic surface and plastic surface. The preconsolidation pressure is the value of yield stress when s=q=0. When C 1 and N(0) are determined, the preconsolidation stress can be determined as follows [3]: In addition, as shown in the later sections, the preconsolidation stress for each test are different and can be changed due to the applied load and wetting processes.
For loading conditions different from triaxial shearing conditions, the objection functions are special cases of Eq. (32) and have simpler expressions, as discussed in the following discussions.

Objective function for elastic zone
Relatively less attention was paid to the calibration of the model parameters in the elastic zone. The elastic parameters   s , and G are considered to be generally of minor importance because of the elastic strains are typically significantly smaller than plastic strain. C 1 in Eq. (24) was not calculated since it does not appear in the original incremental formulation of the BBM. In comparison, calibration of the  and  s for the BBM is nearly the same as the plane-fitting as shown in Fig. 1. The only difference is that specific volume is assumed to be linearly proportional to the net mean stress and suction in the semi-logarithmic scale. The objective function for the elastic behaviour for both isotropic and triaxial stress conditions can be formulated as follows:  Different from the plane-fitting case in which there is a unique C, the C 1 in Eq. (24) under most situations varies for each test and is dependent on the stress history of the soil according to the MSSA. This will be discussed in the Sections 6 and 8 with greater details.

Objective function for elasto-plastic zone
Isotropic loading conditions is a special case of triaxial shearing with q=0. In addition, there is no deviatoric strain under isotropic loading. Consequently, X in Eq.
Consequently, the predicted specific volume v pi can will be calculated using Eq. (25): When a soil is subjected to triaxial shearing conditions, there are both volume change and deviatoric strains. As discussed previously, a nonassociated flow rule was used to calculate the deviatoric strain using Eq. (10), where  is chosen according to Eq. (11). Consequently, all the model parameters in the original BBM can be calibrated using Eq. (32) with the following: and taken as an additional model parameter as suggested by Zhang et al. [15], then the predicted results based upon the calibrated BBM parameter values must be able to match both the volumetric and deviatoric strains simultaneously. Consequently, the object function for the triaxial shearing conditions should be Eq. (29), and where the predicted deviatoric strain ( ) s pi  can be calculated using Eqs. (10) (dilatancy rule) and (13). Under this situation, the deviatoric strain must be calculated using small increments in p, q, and s due to lack of a close form solution for total deviatoric strains. However, the objective functions can be still clearly defined and the components can be easily calculated for the model calibration purposes.

Objective function for K0 loading
The oedometer test is a widely accepted method for studying unsaturated soil behaviour, owing to the simplicity of the experimental setup. In an oedometer test, the lateral stress remains constant at zero throughout the experiment. The condition is a special case of the triaxial loading conditions, which is called, K 0 loading. However, under most cases K0 is unknown. As a result, Hence, the results obtained from oedometer tests are primarily utilized to validate the calibrated model parameters instead of being directly used for model calibration.
Zhang et al. [15] derived an explicit expression to calculate the lateral stress for K0 loading conditions as follows. In the elastic zone, In the elasto-plastic zone, It is noteworthy that the expression for K0 was derived within the critical state framework without any presumptions regarding the soil behaviours. It was formulated based on the zero lateral strain condition that is well satisfied during the oedometer test. Moreover, Zhang et al. [15] developed an optimization approach for objective and straightforward identification of material parameters in the elastoplastic models for unsaturated soils. In the elastoplastic phase, the objective function for Ko-loading is to minimize the errors given by Eq. (29) along the stress path represented by Eq. (36) as follows: It is worth noting that p, q, and in Eq. (36) should exactly follow the stress path as defined by Eq. (36) for K0 conditions. The close-form solution for K0 conditions allows for the direct computation of mean net stress and deviatoric stress from the applied vertical stress and direct use of results from oedometer tests for constitutive modelling purposes. It can simplify the equipment needed for unsaturated soil characterization and significantly reduce the testing time. The time needed for a suction-controlled oedometer test is about 1/10 of the suction-controlled triaxial tests due to much shorter drainage path.
Zhang et al. [15] gave an example to use the results from suction controlled oedometer tests to calibrate the model parameters in the BBM. Zhang et al. [16] used results from constant water content oedometer tests with suction measurements to calibrate the model parameters in the BBM. Each constant water content oedometer test only takes 4-6 hours, which is much less than 1-3 months that is required for the suction controlled triaxial test.
As discussed in Zhang et al. [15], one caution is that the oedometer test is a non-failure test. As a result, the model parameters related to shear strength (M and k in the BBM) using the oedometer test results is only an extrapolation and could be unreliable. It is recommended that other failure tests are used together with the oedometer tests to get a more reliable results for the shear strength. Riad and Zhang [17] combined constant water content oedometer and direct shear tests to calibrate the model parameters in the BBM.

Dealing with Different Stress Histories
To make the test results simple to analyse, most researchers required use of "identical" soil specimens with the same stress histories to perform tests for constitutive modelling purposes [30,32,33]. However, this can be challenging to achieve because unsaturated soil behaviour is notoriously complex, as changes in soil volume and water content can be caused by both stress and suction. Furthermore, unsaturated soil can yield as a result of loading, wetting, drying, or a combination of these processes.
In addition, to develop a constitutive model for simulating unsaturated soil behaviours, numerous laboratory tests are necessary to investigate unsaturated soil behaviour under different stress paths. The resulting data can then be combined to develop a constitutive model capable of predicting soil responses under any arbitrary stress path. In many cases, altering the initial stress conditions is necessary to perform tests with specific stress paths, which has the potential to modify the stress history.
Existing methods for calibrating model parameters for the BBM, as discussed in section 3.1, implicitly assume that all soil specimens have identical stress histories. For instance, in Fig. 2, performing three constant suction isotropic compression tests using specimens with identical stress histories results in the initial yield points H, E, and B falling on the same yield curve, which is used to fit Eq. (3) to obtain p c in the BBM. Achieving identical stress histories for soil specimens requires thoughtful specimen preparation, instrumentation for suction measurements, and lengthy equilibrium time. However, it is only theoretically possible, as loading, drying, and wetting can cause soil to yield and the positions of yield curves to change. For soil specimens with different stress histories, their positions of yield curves will differ, which can lead to misleading results when using existing methods. In Fig. 4, three soil specimens with different stress histories are used for isotropic loading tests under constant suctions s 1 , s 2 , and s 3 . The corresponding v-ln p curves are ABC, D'VF, and G'XI, respectively, and are similar to those in Fig. 2. Using the existing method, curve BVX, composed of three yield points B, V, and X, is considered as the yield curve, but it is incorrect because these yield points belong to three different yield curves LC 1 , LC 2 , and LC 3 , respectively. In contrast, when the three specimens have identical stress histories, the v-ln p curves are ABC, DEF, and GHI, respectively, and have the same yield curve BEH (LC 1 ). However, there are many other possibilities for the relative positions of the yield curves, and the use of existing methods can result in yield curves with significantly different shapes from either BVX or BEH. (c) Three-dimensional representation of volume change of the soil Fig. 4 Test Results for Soil Specimens with Different Stress Histories [28].
When the soil specimens have different stress histories, the test results cannot be used to calibrate the BBM model parameters using existing methods. However, this limitation can be easily overcome by using the MSSA as demonstrated by Zhang et al. [18] and Zhang and Xiao [13]. Fig. 4c shows 3D plot of results in Fig. 4b. It can be seen that AB, D'V, and G'X in Fig.4c do not belong to the same elastic surface and BVX is not a yield curve because these three yield points B, V, and X actually belong to three different yield curves LC 1 , LC 2 , and LC 3 , respectively. However, BC, VF, and XI must fall on the same fixed plastic surface in the v-p-s space as discussed previously, which can be used to best fit Eq. (25) to determine model parameters in the BBM. According to the MSSA, the yield curve is the intercept between the elastic and elasto-plastic surfaces, by letting Eq. (24) equal to Eq. (25), the yield curve can be obtained.
Thu et al. [33] conducted isotropic suctioncontrolled compression tests and constant confining stress drying-wetting tests to determine the BBM LC and SI yield curves using conventional data interpretation methods, as shown in Fig. 5. The MSSA, and Eqs. 24 and (25) were used by Zhang and Lytton [3,4] to reanalysed the Thu et al. [33] data, producing the yield curves shown in Fig. 5. As shown in Fig. 5, the results are vastly different, mainly because the existing methods assume all the soil specimens have the same stress history. However, as discussed in Zhang et al. [28], in order to reach the desired confining stress to run constant confining stress dryingwetting tests, one had changed the stress histories of the soil specimens even if all the soil specimens had the same initial stress histories. Fig. 5 shows that the LC yield curves obtained from the conventional method and the MSSA are vastly different, indicating limitation of the conventional approach.
AB, D'V, and G'X in Fig.4c do not belong to the same elastic surface. When these test results are used to calibrate the model parameters in the elastic zone, different C 1 s should be used as proposed in Eq. (31) in section 5.1.  5. Comparison between yield curves obtained from different methods, from [28] 7 Use of the MSSA to Calibrate the BBM Parameters:

Optimization Strategy
Theoretically, the optimization problem is to find a suitable combination of model parameters to minimize the errors (differences between the measured and predicted results) as defined by the objective function. The optimization can be achieved by a wide range of algorithms [18,21]. The optimization methods can be divided into two groups: (a) direct search methods in which the search strategy is based only on the values of the objective function, similar to the Rosenbrock's [60] and the Simplex methods used by Mattsson et al. [19], and (b) gradient methods, which also require computations of the derivatives of the objective functions. In general, the gradient methods are expected to be more powerful than those using only values of the objective function. On the other hand, the gradient methods are more complex to implement because of the requirements of having an objective function, which always differentiable and which possesses continuous derivatives over the considered domain. The close-rom solutions for the BBM were derived based upon the MSSA as shown in section 4.2. These close-form solutions are second-order differentiable. Zhang and Xiao [13] derived the partial derivatives for the BBM under isotropic conditions. Similar derivatives can be performed easily for triaxial stress states using Eq. (27).
In general, the objective function, F(X), is a highly non-linear function and may possess many local minima. However, the goal of the research is to find a point where F(X) assumes its least value, or the global minimum, inside the considered domain. To achieve the global minimum, a possible strategy could be to run the search algorithm as many times as possible, each time from a different initial condition. Thus, the more runs to cover possible initial conditions, the higher the chance that the global minimum may be found. Another strategy would be to implement a global search algorithm that aims at finding arguments corresponding to the global minimum. In order to find the global minimum of the problem, a combination of quasi-Newton's method and the MSSA objective function (Eqs. [29][30][31][32][33] have been employed in the developed optimization strategy. According to Chong and Zak [34], the Newton method is one of the more successful algorithms for local optimization. If it converges, it has at least a quadratic order of convergence, but if the initial point is not sufficiently close to the solution, the algorithm may not be an adequate one. In this research, a quasi-Newton algorithm, proposed by Ni and Yuan [23], has been used for solving the optimization problem at the local level. The quasi-Newton method combines the advantages of the gradient method with those of the Newton approach. The method uses gradients of the objective functions which yield inverted Hessian matrices, indicating the precision of the obtained parameters. The other advantage of the quasi-Newton method is that there is no need for a resolution of the system of linear functions. To obtain the direction of decrease of the objective function, the quasi-Newton method and gradient method were used to treat inactive and active variables, respectively. Zhang and Xiao [13] provided a detailed description of the optimization algorithm. Besides, comparisons between Newton's method and quasi-Newton's method (BFGS) are also discussed in detail. In a summary, in this research, an objective function that is formulated using the surface approach (MSSA) along with the quasi-Newton's search method is used. Advantageously, the surface approach, which provides an overall best fit to the BBM, is combined with such a powerful search algorithm. The quasi-Newton's method is popular for engineering optimization and programmed into routine office software such as Microsoft Excel® and MatLab®. Consequently, it can be easily accessed and used by ordinary engineers.

Re-analyse the MUSE Benchmark Exercise Data using the MSSA
The procedures described in the previous sections are applied here to demonstrate their simplicity and objectivity for the calibration of model parameters in the BBM. The same experimental data, used in the MUSE benchmark exercise, were reanalysed using the proposed approach to calibrate parameters for the BBM.

Soil properties
The results were taken from the experimental testing program reported in the Ph.D. thesis of Barrera Bucio [24]. The soil samples were collected from natural soils during excavation for the construction of the Rector Gabriel Ferrate Library on the UPC Campus in Barcelona, Spain [2]. Grain size distribution analysis indicated that the soil consisted of 44.5% silt, 39.4% sand, and 16.1% clay which are mainly Illite. The plastic and liquid limits of the soils are 16% and 32%, respectively, and the specific gravity of the soil particles was 2.71. The soil samples were statically compacted by applying all-around pressure of 600 kPa at a water content around 11%, which was approximately 3% dry of optimum. The total suction after compaction was measured as 800 kPa. The samples were then subjected to an initial equalization stage under a matric suction of 800 kPa and low mean net stress values. Additional details about the sample preparation and initial conditions are illustrated in D'Onza et al. [2].

Stress paths and test results
The experimental data used in the MUSE benchmark exercise consisted of six suction-controlled triaxial tests, two isotropic compression tests, and one oedometer test. The matric suction was controlled by means of the axis-translation technique. Fig. 6 presents the stress paths with the identification numbers and the experimental results for isotropic and oedometer tests. Figs. 6 and 7 show the stress paths and experimental results for the suction controlled triaxial experiments, respectively.
Test SAT-1 (Fig. 6a) was an isotropic experiment that consisted of saturating the sample initially by flushing through it (represented by AB in Fig. 6a). Subsequently, the sample was subjected to isotropic loading (BC) until it reached a mean effective stress of 1300 kPa. This was followed by isotropic unloading (CD). Fig. 6d presents the SAT-1 experimental results in e-s-logp space.  Test TISO-1 (Fig. 6b) was a multiple stage suctioncontrolled isotropic test. Initially, the sample underwent isotropic loading (AB) from 0 to 0.6 MPa at a constant suction of 0.8 MPa (AB: s=0.8MPa, p=0→ 0.6 MPa). It was followed by a wetting/drying stage under constant net normal stress of 0.6 MPa (BCD, s= 0.8 MPa →0.01 MPa →0.15 MPa, p=0.6MPa). After that an loading/unloading test stage was performed at a constant suction of 0.15 MPa (DEF: s=0.15MPa, p=0.6 MPa→1.4 MPa→ 0.6 MPa). Finally, a wetting stage at a mean net stress of 0.6 MPa (FG, s=0.15 MPa→0.02 MPa), and an isotropic loading/unloading stage were performed at a constant suction of 0.02 MPa (GHI, s=0.02 MPa, p=0.6 MPa → 2.0 MPa → 0.02 MPa). Fig.  6d presents the TISO-1 experimental results on e-slogp space.
Test EDO-1 (Fig. 6c)  Tests IS-OC-03, IS-NC-06, and IS-NC-12 (Figs. 7a, 7b, and 7c) were suction-controlled triaxial tests in which the mean net stress increased (AB) to 0.3, 0.6 or 1.2 MPa respectively, at a constant suction of 0.8 MPa. The samples were then sheared to failure (BCDE) at constant suction and constant radial net stress, an unloading-reloading cycle was included during shearing for each test. The axial strains versus the applied deviatoric stresses for these tests are presented in Fig. 8.
Test IS-OC-06 (Fig. 7d) was a suction-controlled triaxial test in which isotropic loading and unloading were involved at a constant suction of 0.8 MPa(ABC, s=0.8MPa, p=0 →1.8 MPa →0.6 MPa). The sample was then sheared to failure under the same constant suction and constant radial net stress, with the inclusion of an unload-reload cycle during shearing(CDEF s=0.8MPa, p=0.6 +q/3 MPa, q=0 → 0.8 MPa →0 MPa →failure). The test results are presented in Fig. 8 in Test IWS-OC-01 (Fig. 7e)  Test IWS-NC-02 (Fig. 7f) was a suctioncontrolled triaxial test in which involved  As previously mentioned, a recent benchmark exercise was organized within a "Marie Curie" Research Training Network on "Mechanics of Unsaturated Soils for Engineering" (MUSE) to calibrate parameters for the BBM. In this exercise, experienced constitutive modelers from 7 prestigious teams in unsaturated soils were provided with the same experimental results on an unsaturated soil to calibrate the parameter value for the BBM. The results of this calibration assignment as well as comparisons between the predictions using calibrated parameters from different universities, were published in D' Onza et al. (2015). The calibrated parameters showed high scatter, which indicates the difficulty of calibrating parameters for the BBM.

Calibration of model parameters in the elastic zone
The elastic parameters are calibrated considering results in the elastic zones for all tests. The proposed objective function (Eq. 28) is used for this purpose. Normally there are multiple loading-unloading and wetting-drying cycles in one test. As a result, multiple C 1 values in Eq. (24) might be used for one test to consider the change in stress histories. Using the SAT-01 test as an example, based upon the experimental results, the soil exhibits collapsible soil behaviour. As a result, AB in SAT-01 should be in the elastic zone. from B to C, part of the test results is in the elastic zone and part is in the elasto-plastic zone as demonstrated by the sharp changes in the slope of the test results in Fig. 6d. The yield stress can be approximately determined by the Casagrande method [36] between B and C. From C to D, it is an unloading process, and all the test results should be in a new elastic zone. Consequently, two different C 1 values are needed when the calibration is performed.
The calibrated parameters were =0.0082 and  s =0.0057. These values are in fact close to the values obtained by other universities as shown in table 1. With these values, the experimental results were matched with high accuracy, R-squared of 99.7% is reached. Comparisons between experimental results and predictions are not presented in this section due to its straightforwardness and limited space. These parameters are kept fixed for the following calibration assignments.

Calibration of model parameters in the elastoplastic zone
Zhang and Xiao [13] demonstrated the procedures to calibrate the model parameters in the BBM under isotropic conditions based upon the close-form solution of Eq. (25) according to the MSSA. Zhang et al. [15,16] explain the procedure to calibrate the model parameters in the BBM using suction-controlled and constant water content oedometer test results. In this paper, to simplify the explanation, only the calibration for the triaxial shearing test conditions is discussed. However, all different types of the tests used in the MUSE exercise are used to calibrate the model parameters in the later discussions.
For the triaxial shearing condition, using Test IS-NC-06 as an example, the calibration procedures are as follows: 1. Assume some initial values for ( ) (2016) discussed how to select the initial values for the parameters.
2. Select a yield stress based upon the test results in Fig. 8 for Test IS-NC-06. If the exact initial yield point is uncertain, one can always take a larger value to make sure the selected point is on the virgin surface. For example, at point C when there is unloading, there is irrecoverable axial strain as shown in Fig. 8 give the minimum final sum of all the squares of the difference between the predicted and measured values. Table 3 shows the calibrated parameters from this research (indicated as MST) as well as the parameters calibrated by other universities. The coefficient of determination is also computed. As can be seen in Table 3, the coefficient of determination obtained from this study (97%) is the highest compared with all the other universities. A higher the coefficient of determination represents a better fitting of the experimental results. For example, predictions based upon the parameters provided by UNINN has a R 2 of 93.04%, generating a better fitting of the test results compared with results from other teams. Table 3. Calibrated parameters, standard deviations, and Rsquared from different universities

LC yield curve
As shown in Table 1, all universities in the MUSE exercise provided one preconsolidation pressure * 0 p for as their calibration results. This implies that all of them assumed that the specimens used in the all the tests have the same stress history. Fig. 9 summarized the predicted LC yield curves from all universities. Fig. 9 also shows the results for the initial LC yield curves obtained by the MST team using the method as proposed in . The location of the LC yield curve is obtained by finding the intercept of the movable elastic surface (Eq. 24) and fixed plastic surface (Eq. 25). C 1 in Eq. (24) can be obtained from the initial specific volume of the soil specimen. Once C 1 is known, by letting Eq. (24) equal to Eq. (25) to eliminate the void ratio, the LC yield curve can be found. The preconsolidation pressure can be further calculated by letting s=0. It is found that none of the soil specimens used in the 9 tests have the same initial stress history. It is worth noting that the preconsolidation pressure estimated by the MST team ranges between 0.09 MPa and 0.16 MPa, which was not highly scattered. This indicates that the quality of the tests was relatively well-controlled to have the same stress histories. However, none of them are the same, indicating that it is impossible to implement the conventional approach to obtain the correct initial shape of the LC yield curve. This also explains why the initial yield curves obtained by the 7 universities in the MUSE exercise are significantly different. The conventional approach to determine the LC yield curve has requirements: (1) all the soil specimens have identical stress history, and (2) the yield point can be accurately determined by the Casagrande Method [36]. In practice, none the requirement can be easily met. Accordingly, it is impossible to obtain correct BBM model parameters using the conventional approach as it requires all the soil specimens have the same stress histories. In contrast, the MSSA allows use of "imperfect" results from soil specimen with different stress histories to obtain the correct yield curve. This flexibility can significantly reduce the efforts for producing "identical" soil specimens and simplifies the analysis process.

Comparison between predicted and measured results
The model predictions compared to the experimental results are shown in Figs. 10 through 18. Fig. 10 shows the experimental results and predictions based on calibrated parameters from different universities for suction-controlled triaxial test IS-OC-03 and the corresponding stress path is shown in Fig.7a. The experimental results from the shearing stage are shown in Fig. 10a (in the q- a plane) and Fig. 10b (in the  v - a plane), together with the corresponding model predictions. The shear unloading-reloading cycle was omitted from the model predictions for clarity (this will be done hereafter for all triaxial tests). The dilative behaviour of the soil sample can be clearly seen from Fig. 10b in which the volume of the sample increased as the axial strain increased. The BBM does not consider the soil dilation in its theoretical formulation. Consequently, the model cannot predict this kind of dilative behaviour. Fig. 12a shows predictions from all the seven universities as well as predictions from MST. Obviously, most of the predictions matched the failure strength fairly well. However, most universities under-predicted the axial strain before reaching the failure stage except USTRAT. In general, predictions from MST matched the general trend very well for both pre-failure and at failure stages in Fig. 10a. Fig. 10b indicated that all 7 universities in the MUSE exercise over-predicted the volumetric strain corresponding to axial strain evolution in the early loading stages. Overall, predictions from MST better match the measured behaviour in both Figs. 10a and 10b. Fig. 11 shows the experimental results and predictions based on calibrated parameters from different universities for suction-controlled triaxial test IS-NC-06 and the corresponding stress path is shown in Fig.7b. Only the experimental results from the shearing stage are shown in Fig. 11a Fig. 13a shows all 7 universities in the MUSE exercise underpredicts axial strain, while the MST team slightly over predicts the axial strain, while matching the general trend very well through the whole stress range. Fig. 11b showed that all 7 universities in the MUSE exercise overpredicts volumetric strains, while the MST results can match the test results very well at low axial strain ranges. At high axial strain level, none of the prediction can predict the dilative behaviour of the soil, mainly due to the limitations inherent in the BBM.   Fig. 13 shows the experimental results and predictions based on calibrated parameters from different universities for suction-controlled triaxial test IS-OC-06 and the corresponding stress path is shown in Fig.7d. The soil sample was consolidated under isotropic stress of 1.6 MPa then unloaded to a stress of 0.6 MPa before shearing to failure.
(a) v-lnp curve during the isotropic loading  Figs. 14a and 14b shows the test results under isotropic loading and wetting-drying stages in the v-lnp and v-s planes, respectively. The corresponding model predictions from different groups are also shown. In terms of v-lnp, predictions from most of the universities (except USTRAT) provided reasonable matches to the experimental results. USTRAT team underestimated the yield stress of this soil, which in turn led to over predicting the plastic compression. On the other hand, although DU, ENPC, UNINN, and USTRAT produced predictions matched well with the v-lnp experimental data, they were not able to, simultaneously, provide a satisfactory match for the results on the specific volume changes with suction.
This kind of miss-match can be attributed to a nonappropriate estimation for the LC yield curve shape as discussed previously.
Figs. 14c and 14d show the experimental results from the shearing stage in the a q  − and v a   − planes respectively, as well as the corresponding model predictions.
As shown in Fig. 14c, the MST predictions can match the experimental results at the initial and final shearing stages very well, but largely deviated from the test results before and after the initial yielding. Predictions from the 7 universities in the MUSE exercise overall fit the experimental data better. In terms of va  − , the MST predictions can fit the experimental data very well at the early stage of shearing but fail to predict the dilative behaviour at high deviatoric stress. The MST predictions of va  − relationship for this test is clearly better than those from other research groups.
Experimental results and predictions for the suction-controlled triaxial test IWS-NC-02, using calibrated parameters from various universities, are illustrated in Fig. 15. The corresponding stress path can be seen in Fig. 7f.
The test results for isotropic loading and wettingdrying stages are presented in Figs. 15a and 15b, respectively, depicted in the v-lnp and v-s planes. The corresponding model predictions from various groups are also included. Among all the predictions, the MST team's predictions match the experimental results very well in both the v-lnp and v-s planes, making them the most accurate. Like Figs. 14a and 14b, most other universities predicted the v-lnp relationship well in Fig.  15a, except for USTRAT and DU. Additionally, most universities provided good predictions for the v-s relationship in Fig. 15b, except for USTRAT and ENPC. As a result, only four out of the seven universities were able to reasonably match the experimental results in both the v-lnp and v-s planes. Fig. 15c (in the q- a plane) and Fig. 15d (in the  v - a plane) display the experimental results from the shearing stage, along with corresponding model predictions. However, the experimental results appeared unusual due to the q- a curve's concave shape, which was not observed in any of the other triaxial shearing tests. It is evident from Fig. 15c that all the universities underestimated soil strength and overestimated soil stiffness at the initial loading stages. Despite the curve's concave shape, MST was able to provide an overall good match with the experimental results for all loading stages. In terms of  v - a relationship as shown in Fig. 15d, none of the groups could replicate this dilative response behaviour. However, MST's predictions were the closest to the experimental results for the early stages (i.e., contractive response). plane for stress path BCD in SAT-01 as shown in Fig.  6a. The soil was subject to a loading-unloading cycle under saturated conditions (BCD: s=0MPa, p=0.01MPa→1.3 MPa→ 0.01MPa). Fig. 16 also shows the predicted results for the same stress based on calibrated parameters from different universities as shown in Table 3. Inspection of Fig. 16 showed that the slope of the normal compression line, , as well as the swelling slope,  , were reasonably estimated by all the groups. There were, however, significant variations in the predictions of the yield stress and, consequently, the extension of the elastic zone and the location of the normal compression line. Obviously, both UNITN and ENPC teams overestimated the yielding stress and consequently provided large expansion for the elastic and less plastic volume changes. In contrast, UNINN team underestimated the yield stress value that led to overestimating the volume changes in general. DU, GU, and USTRAT teams, generally, provided proper estimation for the yield stress and better matched the experimental results. The predictions made by the MST is not the best fit but reasonably good. This is consistent with the plane-fitting example, sometimes one must sacrifice the local best fit to achieve a n overall best fit. According to D'Onza et al. [2], none of the teams involved in the MUSE exercise utilized the oedometer test results when selecting parameter values in the BBM. However, in the present study, the oedometer test results (EDO-1) were used to calibrate model parameters in the BBM, based on an explicit formulation proposed by Zhang et al. [15] for calculating lateral stress under K0 conditions. D'Onza et al. [2] noted that predictions by all teams in the MUSE exercise for the oedometer tests typically yielded poorer matches with experimental results compared to isotropic and triaxial tests, likely due to the fact that oedometer test results were not considered in the parameter value determination process. Nonetheless, Fig. 17 demonstrates a strong correlation between the experimental results for EDO-1 and predictions generated using the MST approach, indicating the robustness and mathematical validity of this proposed method for calibrating BBM model parameters, particularly when considering oedometer test results.

Conclusions
The following conclusions can be made from this paper: 1. The development of the BBM was based on the conventional elasto-plasticity approach using incremental formulations. The use of incremental formulation resulted in significant challenges in indepth understanding of the BBM, and subsequent difficulties in model parameter calibration.
2. When the incremental formulation is used, during the calibration process, neither the sequential approach nor the global optimization approach can produce accurate results. This is due to the requirement of an initial condition to calculate the increments necessary for the incremental formulation to initiate the calibration process.
3. The use of incremental formulations also presents difficulties in establishing the objective function required for global optimization, which in turn makes it challenging to achieve accurate results using these methods. The paper demonstrates that even calibrating a simple planar model becomes difficult when using incremental formulations, and the calibration of a more complex model such as the BBM with numerous model parameters and constitutive relations becomes even more challenging. 4. Given the complexity of unsaturated soil behaviour, many tests are typically necessary to calibrate the BBM. Most likely those tests are run with soil specimens with different initial stress histories. However, existing methods often assume that all soil specimens have the same initial stress histories and requires accurate determination of the yield stress. These assumptions, coupled with the use of incremental formulations that presume perfectly correct initial conditions, can lead to significant deviations of the calibrated BBM parameters from the correct results.
5. The MSSA takes full advantage of the uniqueness of the state boundary surface (elastoplastic virgin loading surface), significantly simplifying the process of constitutive modelling. Based upon the MSSA, close-form solutions for the BBM are derived, which makes the establishment of object function for the BBM calibration straightforward. 6. The MSSA can also be used to deal with large amounts of potentially confusing data and synthesize the data into a usable form, which makes the analysis of elastoplastic behaviour of unsaturated soils in a relatively simple way without undue complication. For example, it gives research the flexibility to determine the correct shape of the yield curve and its evolutions during yielding using soil specimens with non-identical stress histories. In addition, based upon the MSSA, constant water content tests and oedometer tests can be easily analysed for constitutive modelling purposes. This significantly reduced the testing time and simplified the equipment needed for unsaturated soil research.