Calibration of a constitutive model for volcanic sands under simple shear conditions

The simple shear response of air-fall volcanic (pyroclastic) soils under both saturated and unsaturated conditions is interpreted through an elastoplastic constitutive model with hydraulic-hardening and porosity-dependent critical state. Extensive experimental data collected under various testing protocols (i.e., simple shear, direct shear, and triaxial tests) enable the theoretical identification of the variability of the fundamental physical properties (i.e., frictional resistance, dilatancy, and water retention behavior) and the identification of a band of admissible values for prescribed confidence levels. In this paper, a constitutive model specifically developed to account for the simple shear loading is adopted and the calibration of the model parameters is performed accounting for such data scatter. The calibration procedure led to a single set of constitutive parameters through 3 main steps. The identification of range of variation for each parameter is carried out using the data derived from laboratory tests. Then, sensitivity analyses are performed on the main physical and constitutive parameters. Mathematical indicators quantifying the difference between the measurements and the predictions are proposed to investigate the role of the material properties and evaluate the model performance. Furthermore, an optimization algorithm is adopted to identify the optimal set of model parameters which best fits all the considered tests (i.e., simple shear tests under saturated and unsaturated conditions and wetting tests). The results imply that the variability of the hydro-mechanical properties must be considered to satisfactorily simulate the constitutive behaviors of the volcanic soils under a variety of simple shear testing regimes.

r Total stress, (kPa) u a Pore air pressure, (kPa) s Matric suction, (kPa) S re Effective degree of saturation, (-) S r Current degree of saturation, (-) S res Residual degree of saturation, (-) h s Saturated volumetric water content, (-) h r Residual volumetric water content, (-) n vg , m vg Shape parameters of the water retention curve, (-) a vg Shape parameter of the water retention curve, (1/kPa) g Stress ratio, (-) g Y Stress ratio at yielding, (-) W State parameter, (-) l g , m g , a g Shape parameters of the dilatancy function, (-) M g The critical stress ratio, (-) N cs InterWcept of the CSL in the n À log r space, (-) k cs Slope of the CSL in the n À log r space, (-) n cs Porosity at critical state, (-) b Model parameter which reflects the effect of the suction on the CSL, (-) e p Plastic normal strain, (-) K Plastic multiplier derived from the consistency condition, (-) D Dilatancy, (-) c p Plastic shear strain, (-) g Y Stress strain at yielding, (-)

Introduction
Shallow landslides are frequently triggered by rainfall in several geo-environmental contexts. A striking example is the 3000 km 2 area of pyroclastic soil deposits in Southern Italy. The response of these natural volcanic air-fall soils along unsaturated steep slope largely depends on their hydro-mechanical coupled behavior, which may also include ''volumetric collapse'' in unsaturated conditions upon wetting, and/or ''static liquefaction'' in saturated condition for undrained shearing. For this reason, experimental and numerical tools are highly needed to identify the material properties responsible for rapid flow sliding [11], as well as estimate the regional landslides susceptibility [12,13,17,25]. At territorial scale, the rainfall-induced landslides are analyzed using different mathematical approaches, but most of the models in the literature adopt the infinite-slope stability analysis to estimate the slope stability [25] and the failure induced by changes in pore water pressure (as Iverson's model, [22], and TRIGRS, [3]. At REV (Representative Element of Volume) scale, the in situ stress-strain conditions, before and at failure, have been investigated through simple shear devices [14,15], which closely approximates the stress-strain paths being the principal (stress/strain) axes free to rotate.
Recent developments extended the use of simple shear testing to partially saturated soils, with the suction that is usually controlled/measured using the axis translation technique and the tests that can be carried out at constant suction or at constant volume. To numerically elaborate the observed soil behaviors, advanced constitutive models reflecting the simple shear conditions while capable to deal with the hydro-mechanical coupling are required. Normally, constitutive models are formulated under triaxial conditions, like the Barcelona Basic model [1,41,43] and the Modified Pastor-Zienkiewicz model [27][28][29]37 . The latter was an extension of the original Pastor-Zienkiewicz model [37] developed for cyclic behavior of saturated soils. To adjust such models to account for the simple shear boundary conditions, 3D generalizations are essential to switch from the space of stress invariants to the manipulations of stress tensor, and furthermore, the stress rotation [36].
In literature some contributions, such as the implementation of the Barcelona Basic Model (BBM) within finite element method, generalize the stress state required to reproduce boundary problems. Gonzalez and Gens [19] showed a coupled flow-deformation finite element analysis of a shallow foundation on an unsaturated loosely compacted silt subjected to water level variations. The behavior of the silt foundation was simulated using the Barcelona Basic Model (BBM), which was implemented in the PLAXIS finite element code. Olivella et al. [36] implemented the BBM in the Code Bright finite element code with different applications in 2D and 3D. Moreover, the Wheeler-Sivakumar model, which derives from the BBM, could be generalized in a similar way of those proposed in literature for the BBM.
However, such generalizations involve intrincated mathematical manipulations and computational efforts. As a simple alternative, constitutive models specifically developed to account for the simple shear loading can effectively avoid such complication while enabling the further understanding of the shearing response of soils [21,38]. In this context, here, a condensed modeling strategy is considered, which intends to reproduce the soil responses under simple shear conditions on the basis of two main premises: (i) the direct use of the control variables that actively control simple shear deformation and failure in the model equations [9,25], and (ii) the explicit incorporation of experimental measurements in the calibration and optimization of the model parameters.
The calibration of the model parameters is carried out through three main steps: (i) identification of parameters' variation range based on experimental results; (ii) sensitivity analysis to identify which parameters have a significant impact on the soil behaviors; and (iii) identification of a single set of parameters through an optimization algorithm to best capture the soil responses under various testing conditions. Hence, special attention is given to the variation of the soil properties and its effect on the model performance by allowing the model parameters to vary within a given range identified from the scatter of the experimental data.

Experimental dataset
The pyroclastic soils derived from the explosive activity of the Somma-Vesuvius volcano (Southern Italy) are generally categorized as sandy silt to silty sand soil (the sandy component ranging between 40 and 60%, the silty component varying from 20 to 60%, and the low gravel component is lower than 22%). They present a high spatial variation of the physical properties [12], but can be clearly subdivided into two main classes (A and B) on the basis of detailed analysis of grain-size distribution and statistical distribution of the physical and mechanical properties [5].
Hereafter, the experimental observations related to remolded specimens are detailed. As well reported by Bilotta et al. [5], soil strength features have been comprehensively investigated through triaxial tests and direct shear tests in saturated and unsaturated conditions. Recent developments extended the potential of simple shear tests to partially saturated soils, using the axis translation technique to control/measure the suction of the specimens. Simple shear tests overcome some limitations of both the triaxial apparatus and direct shear device [8,39] as the rotation of the stress-strain principal axes can be reproduced. The tested specimens were firstly consolidated under K0 conditions. Then, for the saturated samples, the simple shear tests (i.e., SSP0115, SSP0215, and SSP0315 in Table 1) were carried out under drained condition by applying a shear rate of 0.002 mm/min [14]. The simple shear tests on unsaturated specimens (i.e., SSRPSF03a, SSRPSF03b, SSRPSG23, and SSRPSG24 in Table 1) were performed using the same procedure at constant suction and constant normal stress [31]. Generally, the saturated specimens experienced the contractive behavior coupled with a soil hardening behavior, while the unsaturated specimens showed the transition from contractive to dilative response associated with a slightly hardening behavior [2].
Simple shear equipment was also adopted to investigate the mechanical characteristics of the soil upon wetting. The wetting tests (i.e., SSRPSG03 and SSRPSG14 in Table 1) were performed at different initial stress levels and the suction is applied at either the bottom or the top of the specimens using several suction variation rates and measured on the other side. The stress ratio (i.e., shear stress divided by effective vertical stress at failure) versus the suction measured at failure was obtained for different tests [33], and for a suction variation rate equal to -0.1 kPa/h, the failure occurred at a stress ratio between 0.80 and 1.10.

Constitutive model
The behavior of pyroclastic soils in simple shear condition under saturated or unsaturated regimes and upon the transition from unsaturated to saturated condition was analyzed through an elastoplastic strain-hardening constitutive model proposed by Buscarnera and di Prisco [10] with specific adjustments made to the dilatancy law for the studied pyroclastic soils under the simple shear stress conditions. The Bishop effective stress [7], specifically adjusted to the constitutive modeling of geomaterials [23,40], is adopted here: where r is the total stress, u a is the pore air pressure, s is the suction, and S re is the effective degree of saturation defined as S re ¼ S r À S res =1:0 À S res where S r is the current degree of saturation, and S res is the residual degree of saturation. The effective degree of saturation is estimated using the van Genuchten retention model [42], expressed in terms of volumetric water content (i.e., S r ¼ h=n): where h s is the saturated volumetric water content, h r is residual volumetric water content, and a vg , n vg and m vg ¼ 1 À 1=n vg are shape parameters of the water retention curve.
A frictional yielding criterion is adopted by assuming a negligible cohesive response [25], formulated as: where g Y is the stress ratio at yielding, s is the shear stress, and r 0 is the effective normal stress.
A porosity-dependent non-associated flow rule is adopted by enhancing the dilatancy function proposed by [24] with the concept of state parameter w [4,26], expressed as: where e and c are the normal and shear strain with the superscript p representing their plastic parts, l g , m g , and a g are shape parameters of the dilatancy function, M g is the critical stress ratio governed by the friction angle u (i.e., Fig. 1 Schematics of the suction-dependent critical state conditions and consequent changes in dilatancy M g ¼ tanu), and the state variable, w, quantifies the distance between the current porosity and the porosity at the critical state, defined as: where the porosity at the critical state under the saturated regime is expressed as: in which N cs and k cs represent the intercept and slope of the CSL in the n À log r 0 space. It is worth noting that the incorporation of the state variable in the dilatancy equation enables the consideration of porosity-dependent features of sandy soils. Specifically, for a loose soil with a positive w (i.e., Point A in Fig. 1), a contractive response is predicted (i.e., D [ 0), while for a dense soil with a negative w (i.e., Point B in Fig. 1), a dilative response is predicted (i.e., D\0). The model has the capacity to capture the transition from contractive (i.e., D [ 0) to dilative response (i.e., D\0). Furthermore, as suggested by the experimental evidence [18], the presence of suction could cause an upward shift of the CSL leading a contractive response under saturated conditions to dilate while under unsaturated conditions. Like Point A in Fig. 1, w takes positive values under saturated states but it becomes negative after drying. In this context, the critical state line needs to be adjusted to incorporate the effect of suction. Although the work by Gallipoli et al. [18] suggests that the critical porosity varies exponentially with suction, a linear relationship is simply adopted to account for the impact of suction on the critical state line, being: where b is a fitting parameter. And the plastic strains can be computed as: where K is the plastic multiplier derived from the consistency condition. It is worth noting that since the simple Trial-and-error procedure -k -Trial-and-error procedure --CSL N cs -2S S sat 2SS unsat n cs ¼ N cs À k cs ln r n cs À ln r k cs -Shear strength / 0°3 SS sat 4 SS unsat 11 DS sat 3 DS unsat 7 TX sat 5 SS wet shear stress ratio quantitatively reflects the degree of stress rotation [34,35], the proposed dilatancy formula reflects, only in an implicit manner, the effect of principal stress rotation on the development of the plastic strain. It must therefore be noted that full consideration of possible anisotropic characteristics of the plastic flow emerging from principal stress rotation would require a more advanced treatment accounting for multi-axial stress effects and the possible non-coaxiality of the plastic flow rule [20,21].
Shear-hardening is captured through the following strain-hardening law [25]: where k is a hardening constant. As the material response approaching the critical state, g Y gradually coincides with M g , which effectively prevents the growth of yield surface at the critical state. Finally, the elastic response is modeled through pressure-dependent hypoelasticity with the Young's modulus E and shear modulus G expressed as: where n E and n G are constant power law coefficients; E r and G r are the values of the elastic moduli at the reference stress r r (normally being 1 kPa).

Identification of parameters' variation range
The variation range of the model parameters was accomplished based on a series of laboratory tests performed in previous studies.   Based on the statistical distribution of the experimental data, five confidence intervals were identified to comply the variation of soil properties associated with SWRC, dilatancy, and shear strength. Areas bounded by 10% 25%, 50%, 75%, and 95% confidence levels are analyzed which suggests the probability of a new observation lying within a predicted interval. The procedure adopted for the identification of parameters' variation range under a given confidence level is detailed below and Table 2 reports the considered equation, the minimum number of experimental tests, and the variable plane in which the ranges were calibrated.
The SWRC was described through three parameters, which were obtained using two wetting tests in simple shear condition. The ranges of variation of the SWRC parameters were identified considering only the experimental stage in which the shear stress (34 kPa, or 51 kPa) and net vertical stress (30 kPa, or 50 kPa, respectively) were maintained constant, while the suction was lowered at a suction rate of -0.1 kPa/h (see Sect. 2). In this context, the parameters in Eq. 2 were determined and are listed in Table 3 with the boundary values describing the various confidence intervals, as shown in Fig. 2. By analyzing the data, it was possible to define a single value of a (i.e., a single Air Entry value, AEV, which is the suction value that corresponds to the entry of air in the largest pores of the soil) and a narrow interval of variability for the parameter n.
The ranges of variation of the elastic parameters (E r , G r , n G and n E ) were estimated using triaxial tests in drained condition on saturated specimens. The parameters E r , G r , n G and n E were obtained by fitting Eq. 10 with the experimental data obtained from the triaxial tests.
The critical state parameters N cs and k cs were constrained using the n À log r 0 plane. Since none of the TX, DS and SS tests fully reached the critical state, the precise  determination of the location of the Critical State Line (CSL) was impossible. Therefore, the CSL was approximately evaluated by hypothesizing that the Critical State is reached at large deformations (e.g., for TX tests it was supposed that CSL was reached at axial strain equals to 40%). For the SS tests, the CSL was defined by extrapolating the experimental results for an extra 5% of shear strain. Here, according to previous works [30], the parameters at critical state are k cs ¼ 0:031 and N cs ¼ 0:78. However, the confidence intervals could not be properly identified due to the uncertainty about the critical state. On the other hand, the parameter M g , being a function of friction angle /, was calibrated on the basis of the results obtained in TX, DS and SS tests. For each test reported in Table 1, by assuming negligible cohesion, the friction angle can be computed as tan/ ¼ s=r 0 using the stress values at failure. As elaborated by Bilotta et al. [5], the saturated conditions would cause substantial differences in the measured friction angle. Therefore, it would be essential to consider the variation of friction angle in the numerical assessments of material responses. The highest friction angle was observed in simple shear tests under unsaturated regimes; conversely, the lowest friction angle was observed for unsaturated specimens in direct shear condition.
To incorporate the suction impacts on the shear strength, the variability of friction angle was evaluated in relation to the normal effective stress. An exponential correlation was identified in the / À r 0 plane (Fig. 3), and the confidence levels specified above were determined correspondingly to account for the variability (Table 3).
where / 0 , c 1 and c 2 are model constants.
Although there is evidence of nonlinear trends in the failure envelope, maybe related to nonlinear water retention effects on the apparent cohesion, capturing them has not been considered a priority in this study.
The variation of the friction angle also influenced the calibration of the dilatancy parameters. They were calibrated through two simple shear tests, one saturated (i.e., SSP0115) and one unsaturated (i.e., SSRPSG24). The parametersl g , a g and M Ã g were evaluated by fitting the measured dilatancy with Eq. 4, while m g was obtained through the reversal formula of Eq. 4, M g , where M g was estimated using the fitted parameter related to shear strength. The estimated range of l g , a g and m g are shown in Table 3, and the confidence levels of dilatancy are also plotted in Fig. 4. Finally, the hardening parameter (k) was estimiated through a trial-and-error procedure, being 0.06, and then further adjusted using the sensitivity analysis and optimization procedure described in the next sections.

Sensitivity analysis
In the previous section, it was outlined how the variability in soil properties influences the choice of the optimal constitutive parameters. Here, a local sensitivity analysis was performed by turning the value of the studied parameter, while keeping the remaining ones constant. Local sensitivity analysis is performed to assess the local impact of material variations on the model response within the captured fluctuating range. This method is particularly suitable when the importance of constitutive features on the material response is not clear, even if it requires computational iterations on each individual parameter over a defined range to clarify its significance on the systematic behavior.
To provide an objective and quantitative evaluation of the model performance, an error function was defined to quantify the difference between measurements and modeling results. Given n unknown variables named as F i with i = 1,.., n, stress or strain, the error function Err F i was defined as the relative distance between the measured value, F exp i , and the modeled value, F mod i , computed as: The error function was evaluated at 10 points with specified values of shear strain (cÞ for the simple shear tests. Specifically, the shear strain values were selected with a step of 1% in the range 0\c\5% and with a step of 5% for c [ 5%. The points considered are marked in the following plots as solid red circles. It was expected that the error function may vary along one single test, and differently from one test to another, depending on the type of test, confining stress, and specimen condition. The sensitivity analysis was performed on 4 simple shear tests (2 saturated and 2 unsaturated), and 2 wetting tests in simple shear condition. The model performance was evaluated in the s À c and e À c planes, and for the unsaturated tests, the S r À c plane was also considered.
The mean error in shear stress Err s m;j and the mean error in normal strain Err e m;j were evaluated for each test j ð Þ, and the mean error in the degree of saturation Err S r m;j was also evaluated for the simple shear unsaturated tests and for the wetting tests. Moreover, for each test, the mean total error (Err tot,m ) was estimated as the sum of the mean error in the shear stress, the normal strain, and the degree of saturation: The first group of parameters analyzed was related to the friction angle, which strongly influences the model performance. Simulations of the two simple shear tests on saturated specimens were performed for all the confidence levels identified in Sect. 4 using the data reported in Table 3. As shown in Fig. 5, the friction angle influences the model performance either directly through M g (as indicated in the s À c plane) or indirectly through dilatancy (as indicated in the e À c plane). The variability in the friction angle induces significant changes in terms of maximum shear stress and shear strain reached at large deformation. The slope of the s À c curve at low shear strain (c\5%) was not heavily affected by the variability in friction angle (Fig. 5a and e).
The test SSP0115 exhibits the lowest Err s m (being 8%) at the lower bound of the 95% confidence level, while the lowest Err e m (being 10%) corresponds to the upper bound of the 50% confidence level. On the contrary, the test SSP0315 exhibits the lowest Err s m (being10%) at the two lower bounds of the confidence level equal to 75% and 50%, but the lowest Err e m (7%) is observed at the lower bounds of the confidence level equal to 75%. The range with the lowest Err tot; m identified by simple shear tests in Similarly, as shown in Fig. 6, the effect of friction angle on unsaturated specimens was also significant in the s À c and e À c planes, while the error Err S r is almost constant since it is basically controlled by the shape of the SWRC. The observations done for the saturated tests are still valid here, and the friction angle value mainly influences the reaching of maximum shear stress and the maximum normal strain at larger deformation. The volumetric behavior was well captured by the majority of the values attributed to the parameters related to the shear strength. However, the error related to the normal strain was very high Þat shear strain between 25 and 35% for the SSRPSG24 test (Fig. 6e). This is due to the small absolute values used in computing the error. In the SSRPSG23 test, the error related to the normal strain is always lower than that observed for the SSRPSG24 test but associated with an increasing trend with shear strain (Fig. 6). The range with the lowest Err tot; m identified by the simple shear unsaturated tests was located between the curve fit and the upper bound of the confidence levels 25%.
A sensitivity analysis for the parameters related to dilatancy law (a g , l g and m g ) was then performed. Under saturated conditions, those parameters only influence the model prediction in the À c plane (Fig. 7b, d). Although the mean total errors for tests SSP0115 and SSP0315 were included in a narrow range for the various considered confident intervals, the best model performances (i.e., those corresponded to the lower mean total error) were obtained using the set of parameters corresponding to lower limits of confidence level 75%. The best model predictions for saturated samples were obtained using the parameters reported in Table 4, which corresponds to a narrow range around the lower bound of the confidence level 75%. Similarly, dilatancy parameters mainly influenced the model predictions in the e À c plane also in unsaturated condition. For each test, the bounds of confidence level with low mean total errors were identified and the best model predictions were ensured by the parameters corresponding to lower bounds of the confidence levels 75% and 95%. The lowest mean total error corresponds to the parameters related to the lower bound of the confidence level 75% (Table 4, in the column for unsaturated shearing).
A sensitivity analysis on the hardening parameter was performed with k varying between 0:01 and0:11, which was suggested as reasonable by the trial-and-error procedure. The highest Err s m . and Err e m were observed at low shear strain (c 10%), because for a decrease in the hardening parameter the maximum curvature point shifted toward smaller shear strain (Fig. 8a and e).
Thus, smaller values of k led to an increase in the slope of the elastic part of s À c curve. The maximum shear stress was independent on the k value, contrary to the maximum normal strain, i.e., the larger k, the higher the maximum normal strain. The best fit between simulations and measurements was reached at k ¼ 0:07. As already noticed in saturated condition, the paramater k influenced the model prediction in both the s À c and e À c planes. The best model performances were obtained at k ¼ 0:05 for SSRPSG23 and k ¼ 0:07 for SSRPSG24. The value of k, which provided the minimum total error for the two tests, was 0:07.
The sensitivity analysis for the parameter b was carried out only with reference to the two simple shear tests under unsaturated conditions. The same procedure used for k was also used for the b parameter and the best fitting performance was obtained at b ¼ 0:02. The parameter b mainly influenced the model performances ine À c, while only negligible influences were observed in s À c plane (Fig. 9). The highest Err s m and Err e m were observed at larger shear strain.
Since the measured SWRC responses were well captured by the linear-fitting curve and the variation in the shape of SWRC has very limited impacts on shearing responses, the variability of SWRC was not considered in the sensitivity analysis and the values corresponding to linear-fitting curve were adopted.
The effect of soil properties variation was also a critical issue simulating the behavior of pyroclastic soils upon wetting; thus, the sensitivity analysis was also performed on the two wetting tests reported in Moscariello [32] (i.e., SSRPSG03 and SSRPSG 14). The wetting tests were performed under constant vertical stress and shear stress under a constant suction rate of -0.1 kPa/h [16]. The sensitivity analysis was carried out for the parameters related to shear strength, dilatancy law, and hardening parameter k, using the same ranges of variation identified by the shearing tests on saturated and unsaturated specimens. Figure 10 summarizes the mean prediction errors of s and e at each confidence level of the studied parameters for shearing tests under saturated and unsaturated conditions and wetting tests. The parameters' range of variation with the lowest mean total error were identified within the marked gray zone, corresponding to the minimum total error (the sum of Err s and Err e ). Comparing the three testing regimes (i.e., saturated shearing, unsaturated  different considering saturated, unsaturated and wetting tests, which are consistent with the experiment observations [5]. The parameters with the best model performances are summarized in Table 4.

Optimization
An optimization procedure was applied to improve the model performance on four simple shear tests (2 SS sat and 2SS unsat ). The sensitivity analysis outlined that shear strength and plasticity parameters induced high changes in model behavior. Thus, the analyses of the tests were done by changing, under a controlled range of variation identified in Sect. 2, the parameters related to shear strength and plasticity. Specially, five parameters were optimized: / 0 , c 1 , c 2 (shear strength); m g and k (plasticity). These parameters correspond to those with a different optimal value for each type of test considered.
The optimization was carried out through the Optimization Toolbox of MATLAB, which uses an iterative procedure, searching a set of input parameters that minimizes the difference between the model outcome and the experimental data. These differences can be called as error function or objective function, here defined as follows: where Err tot;m;j is the mean total error of the test j as defined in Eq. 13. The Optimization Toolbox of MATLAB contain Multiobjective minimizers, which attempt to either minimize the maximum value of a set of functions, or to find a location where a collection of functions is below some The description of the specific features of the solver function and the differences among those available in the Optimization Toolbox of MATLAB or in other similar codes are out the purposes of the paper; thus here, only the main features of the solver and the constraints adopted are described. Some algorithms (i.e., fminsearch, lsqnonneg, fmincon), which find the optimum solution of multivariate and nonlinear functions, were tested, and some of these provided good results. However, the solver, supplying the best performances, was ''fmincon,'' which finds minimum of constrained multivariable function using derivative-free method. ''fmincon'' is a gradient-based method that is designed to work on problems, where the objective and constraint functions are both continuous and have continuous first derivatives. The algorithm used inside the solver function is named ''interior points'', which satisfies bounds at all iterations and can handle large, sparse problems, as well as small dense problems. This algorithm was selected because needs of low memory usage and exhibits the ability to solve large problems quickly. Moreover, it is possible imposing the parameters ranges for the confidence levels identified in Sect. 4. The bound constraints of / 0 , c 1 , c 2 and m g were chosen to be the bounds corresponding to the 95% confidence level, which are shown in Table 3. The bound constraints of k were 0.01 and 0.09. The values of / 0 , c 1 , c 2 , m g and k, resulting from the optimization procedure, are shown in Table 5.

Model performances
Both the simple shear saturated and unsaturated tests, and the simple shear wetting tests were simulated using the unified optimized parameters reported in Table 5.
The model parameters related to friction angle obtained from the optimization procedure are comprised between the lower bounds of 25% and 75% confidence level. The value of the hardening parameter is equal to 0.06, which was already identified as the optimum value to predict the behavior of the simple shear saturated tests. The optimum value of m g corresponded to the 50% confidence level.
The general agreement between the experimental results and the constitutive modelling is discussed with reference to s À c and e À c planes for simple shear saturated and unsaturated tests, and with reference to c À s and e À s planes for the simple shear wetting tests.
The model prediction well reproduces the trends of all the experimental results. Interestingly, the numerical results well reproduced the pyroclastic behavior in the s À c plane, and very low differences were found in the e À c plane for simple shear tests, c À s and e À s planes for the simple shear wetting tests. The best numerical simulations in both planes were obtained for SSP0315 (saturated), SSRPSG24 (unsaturated) and SSRPSG14 (wetting test). Low values of Err e were observed at shear strain higher than 10% (Fig. 11b, d and f). In the s À c plane, the Err s is lower than 12% for SSP0115 (being 10%) and SSRPSG24 (being 11%) tests, while is lower than 45% for the SSP0315 (being 42%) and SSRPSG23 (being 38%) tests. The median error related to the degree of saturation was in all tests lower than 10% (Err S r m =7%, in SSRPSG23 and Err S r m =0.12% in SSRPSG24).

Concluding remarks
The paper shows a strategy to calibrate a constitutive model using a large database of experimental tests, and taking into account the variability of soil physical and mechanical properties. The soil tested is a pyroclastic soil originated from eruptive activities of the Somma-Vesuvius volcano. The model adopted is specified to reprod uce simple shear condition and several saturation condition, i.e., saturated, unsaturated condition and the transition from unsaturated to saturated condition. The model was calibrated through 28 laboratory tests. The laboratory tests included tests performed on saturated and unsaturated specimens and four mechanisms of failure were considered, i.e., triaxial, simple, direct shear and wetting in simple shear conditions. This large amount of tests allows well assessing the parameters' range of variations due to the large variation of soil properties. Based on the variability of the soil properties, five different confidence levels (and their boundary values) were identified for each variable. The confidence levels were estimated for the parameters related to shear strength and volumetric behavior (i.e., dilatancy law), while a similar analysis was impossible for the parameters related to critical state. For each confidence level, the parameters related to its lower and upper bounds were identified.
A sensitivity analysis was carried out to investigate the impact of the changes in friction angle, dilatancy and hardening parameters on model performances, which were evaluated using an error function adjusted from previous research works. The analysis also allowed the identification of the mechanical properties of the pyroclastic soil that most impacted the performance of the constitutive model. Fig. 11 Comparison between measurements and numerical results obtained using the calibrated parameter set proposed in Table 3: a, b shearing tests under saturated conditions; c, d shearing tests under unsaturated conditions; e, f wetting tests; g prediction errors The model performances were estimated by simulating two simple shear tests on saturated specimens, two simple shear tests on unsaturated specimens and three wetting tests in simple shear condition. The bound or the bounds of the confidence levels, which guaranteed the minimum total mean error, were identified and later used as the constraint during the optimization procedure.
The sensitivity analysis carried out for the simple shear tests outlined that the parameters, which largely influenced the model performance under saturated and unsaturated condition, were those related to shear strength, i.e., the friction angle. The choice of friction angle parameters influenced the model performances in terms of magnitude of the maximum shear stress and normal strain both in saturated and unsaturated condition. The mean total errors were higher than the ones observed by varying the dilatancy law parameters. The hardening parameter also conditioned the model performances in s À c and e À c planes, but increasing k value, the point of maximum curvature of the s À c curve was shifted toward smaller shear strain values. This parameter played a significant role in the wetting tests in simple shear condition, and this can be attributed to the low stress strain and stress shear reached upon the shearing stage.
The optimization procedure, accomplished using the Optimization Toolbox of MATLAB, allows a single set of constitutive parameters to be obtained, which provided satisfactory simulation of the hydro-mechanical soil response in simple shear condition under several saturation condition and failure mechanisms. The model was calibrated and validated for one soil in three different soil conditions. Further investigations could extend the application of the calibration and optimization procedure to other soils.
Funding Open access funding provided by Università degli Studi di Salerno within the CRUI-CARE Agreement.
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/.