Jie Yang

Dr Yin-Fu Jin

Dr Zhen-Yu Yin
email: zhenyu.yin@polyu.edu.hk

Farid Laouafa

Pierre-Yves Hicher

Keywords: Granular soils, internal erosion, hydro-mechanical coupling, parameter identification, optimization

Due to the complexity of the hydro-mechanical behavior of soils subjected to internal erosion, a high number of parameters is usually required for the erosion models and the constitutive models. This aspect makes it difficult to determine by trial-error their relevant values from laboratory tests. To address this issue, an efficient optimization-based procedure for identifying the parameters of a recently developed hydro-mechanical model for internal erosion using an enhanced backtracking search algorithm (so-called MBSA-LS) has been proposed. The MBSA-LS incorporates two points: (1) modifying the mutation of the original BSA (Backtracking Search Algorithm) and (2) incorporating an efficient differential evolution (DE) as a local search to improve the optimization performance. A mono-objective framework with six different criteria has been proposed to identify the parameters related to the interlocking effect and the erosion process. The proposed procedure was successfully applied to identify the parameters from the erosion tests of Hong Kong-Completely Decomposed Granite mixture (HK-CDG). All results demonstrated that coupling the MBSA-LS and the hydro-mechanical erosion model could efficiently solve the issue of parameter identification accounting for both the mechanical behavior and internal erosion.

/ 29

Identifying the parameters of a hydro-mechanical model for internal erosion occurring in granular soils by using an enhanced backtracking search algorithm

Introduction

Internal erosion is a significant issue in civil and environmental engineering, which occurs when fine particles are plucked off by hydraulic forces and transported through the coarse matrix by seepage flow [START_REF] Bonelli | On the modelling of internal soil erosion[END_REF][START_REF] Fell | Time for development of internal erosion and piping in embankment dams[END_REF][START_REF] Wan | Investigation of rate of erosion of soils in embankment dams[END_REF]). Due to the loss of a part of the solid phase, the mechanical properties of the soil are progressively modified and degraded. Statistical analyses of accidents in embankment dams indicate that one of the major causes of accidents in embankment dams is internal erosion [START_REF] Fell | Time for development of internal erosion and piping in embankment dams[END_REF][START_REF] Foster | The statistics of embankment dam failures and accidents[END_REF][START_REF] Muir Wood | The magic of sands-the 20th Bjerrum Lecture presented in Oslo, 25 November 2005[END_REF][START_REF] Rönnqvist | On the use of empirical methods for assessment of filters in embankment dams[END_REF][START_REF] Xu | Breaching parameters for earth and rockfill dams[END_REF][START_REF] Zhang | Seepage failure mechanism of the Gouhou rockfill dam during reservoir water infiltration[END_REF]; Zhang et al. 2009a). Furthermore, recent studies have shown that internal erosion is also an important issue in geotechnical structures, such as land subsidence due to water piping induced erosion [START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF], lateral displacement induced by erosion during jet grouting [START_REF] Shen | Estimation of lateral displacement induced by jet grouting in clayey soils[END_REF], ground surface settlement induced by erosion due to tunnel leakage [START_REF] Wu | Identification of tunnel settlement caused by land subsidence in soft deposit of Shanghai[END_REF], and landslides or slope instability induced by fines migration during rainfall [START_REF] Crosta | On slope instability induced by seepage erosion[END_REF][START_REF] Hu | Seismic precursor to instability induced by internal erosion in loose granular slopes[END_REF][START_REF] Lei | Numerical investigation of rainfall-induced fines migration and its influences on slope stability[END_REF]. Over the last decades, internal erosion has been widely studied by laboratory tests using primarily the conventional rigid-walled permeameter [START_REF] Kenney | Internal stability of granular filters[END_REF][START_REF] Rochim | Effects of Hydraulic Loading History on Suffusion Susceptibility of Cohesionless Soils[END_REF]Sibille et al. 2015a;[START_REF] Skempton | Experiments on piping in sandy gravels[END_REF][START_REF] Sterpi | Effects of the erosion and transport of fine particles due to seepage flow[END_REF]; [START_REF] Wan | Assessing the potential of internal instability and suffusion in embankment dams and their foundations[END_REF]. The eroded mass, the hydraulic gradient and the flow rate were monitored to determine the likeliness of seepage induced instability by incorporating the post-erosion analysis of the particle size distribution.

Recently, more attention has been paid to the coupled hydro-mechanical behavior during the erosion process; in particular the role of the stress ratio under which the erosion process was active has been emphasised [START_REF] Bendahmane | Experimental parametric study of suffusion and backward erosion[END_REF][START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF]Ke and Takahashi 2014b;[START_REF] Marot | Assessing the susceptibility of gap-graded soils to internal erosion: proposition of a new experimental methodology[END_REF][START_REF] Slangen | A Flexible Wall Permeameter for Investigating Suffusion and Suffosion[END_REF].

Based on experimental findings, numerical analyses at the scale of the representative elementary volume have been proposed to describe the erosion process, either within the framework of continuum mechanics by the finite element or the finite difference methods, or the material point method [START_REF] Cividini | Finite-element approach to the erosion and transport of fine particles in granular soils[END_REF][START_REF] Papamichos | Volumetric sand production model and experiment[END_REF][START_REF] Stavropoulou | Coupled wellbore erosion and stability analysis[END_REF][START_REF] Vardoulakis | Hydro-mechanical aspects of the sand production problem[END_REF]Yang et al. 2019aYang et al. , 2019b;;[START_REF] Yerro | Modelling internal erosion with the material point method[END_REF][START_REF] Zhang | A thermodynamics-based model on the internal erosion of earth structures[END_REF], or within a discrete framework by coupling the discrete element method with the computational fluid dynamics [START_REF] Lominé | Modeling of fluid-solid interaction in granular media with coupled lattice Boltzmann/discrete element methods: application to piping erosion[END_REF][START_REF] Mansouri | Numerical simulation of the quicksand phenomenon by a 3D coupled Discrete Element -Lattice Boltzmann hydromechanical model[END_REF][START_REF] Reboul | A statistical analysis of void size distribution in a simulated narrowly graded packing of spheres[END_REF][START_REF] Sari | Investigation of internal erosion processes using a coupled dem-fluid method[END_REF]Sibille et al. 2015b;[START_REF] Zhang | Three-dimensional DEM modeling of the stress-strain behavior for the gap-graded soils subjected to internal erosion[END_REF][START_REF] Zhao | Coupled CFD-DEM simulation of fluid-particle interaction in geomechanics[END_REF]. Since the continuous approach permits numerical simulations to be done at the entire engineering structure scale, it is undoubtedly preferable to operate within this framework and to improve the mechanical response of the soil by monitoring the porosity and the evolution of the fines content throughout the process of internal erosion. However, modelling the hydro-mechanical coupling induced by internal erosion is still quite rare at this scale. More recently, Yang et al. (2019c) proposed a numerical model within the framework of continuum mechanics for investigating the coupled hydro-mechanical behavior of soils subjected to internal erosion.

The approach was able to describe the process of internal erosion, the deformations of the specimen during erosion, and the degradation of the mechanical properties of the eroded specimen. The proposed numerical model was applied to assess how internal erosion impacts the safety of earthen structures (Yang et al. 2019d). However, adding refinement and accuracy to model results in the introduction of complexities along with additional parameters. This requires determining the parameters before the model is applied to actual engineering scale problems. Since the parameters are not easy to quantify with conventional laboratory tests, an efficient procedure for the parameter identification should be developed.

Recently, inverse analysis methods by optimisation have been successfully used in the geotechnical field (Jin et al. 2016a;[START_REF] Jin | Identifying creep and destructuration related soil parameters by optimization methods[END_REF][START_REF] Jin | Identifying parameters of easily crushable sand and application to offshore pile driving[END_REF]Jin et al. 2019a;Jin et al. 2019b;Jin et al. 2019c;[START_REF] Levasseur | Soil parameter identification using a genetic algorithm[END_REF]Yin et al. 2018b;[START_REF] Yin | Identifying controlling soil delayed behaviour from laboratory and in situ pressuremeter testing[END_REF]; apparently these methods produce a relatively objective determination of the parameters for a given soil model, even for those that have more of a mathematical than a physical meaning.

The existing optimization techniques can be divided into two main categories: (1) deterministic optimization techniques, such as gradient-based and simplex algorithms [START_REF] Calvello | Selecting parameters to optimize in model calibration by inverse analysis[END_REF][START_REF] Papon | Single-and multi-objective genetic algorithm optimization for identifying soil parameters[END_REF][START_REF] Yin | Identifying controlling soil delayed behaviour from laboratory and in situ pressuremeter testing[END_REF], which deal with a single solution and are focused on reaching local minima since they start the search procedure with a first trial value (often chosen randomly in a specific space). If this supposed value is not close enough to the global minimum solution, the solution is likely to be trapped in a local minimum.

(2) Stochastic optimisation techniques, such as evolutionary algorithms [START_REF] Levasseur | Soil parameter identification using a genetic algorithm[END_REF][START_REF] Pal | Calibration of a constitutive model using genetic algorithms[END_REF][START_REF] Papon | Single-and multi-objective genetic algorithm optimization for identifying soil parameters[END_REF], Particle Swarm Optimiser (PSO) algorithms [START_REF] Knabe | Identification of constitutive parameters of soil using an optimization strategy and statistical analysis[END_REF]Zhang et al. 2009b) and differential evolution algorithm [START_REF] Zhao | Inverse analysis of deep excavation using differential evolution algorithm[END_REF], the artificial bee colony algorithm (ABC) [START_REF] Kang | An artificial bee colony algorithm for locating the critical slip surface in slope stability analysis[END_REF][START_REF] Kang | Artificial bee colony algorithm with local search for numerical optimization[END_REF], require high computational resources. Among these optimisation algorithms, the stochastic algorithm has attracted more attention given its greater capacity to solve complex problems in the geotechnical field (Jin et al. 2016a;Jin et al. 2016b;Jin et al. 2016c;Jin et al. 2016d;[START_REF] Levasseur | Soil parameter identification using a genetic algorithm[END_REF][START_REF] Papon | Single-and multi-objective genetic algorithm optimization for identifying soil parameters[END_REF]Yin et al. 2016b;Yin et al. 2016c). Among numerous stochastic optimizations, the Backtracking Search Algorithm (BSA) has a strong ability of exploitation and exploration due to its mutation and crossover strategies in which the historical population information is used to guide the search. However, the BSA has not yet been applied to identify the parameters of a coupled hydro-mechanical model of internal erosion.

Therefore, this study proposes an efficient optimization procedure for identifying the parameters of a hydro-mechanical model of internal erosion for granular soils. For this purpose, a recently proposed coupled hydro-mechanical model of internal erosion for granular soils has been adopted for the simulations. A novel backtracking search algorithm enhanced by differential evolution (DE) as a local search has been used to find the optimal solution. A mono-objective parameter identification procedure that combines the hydromechanical model of internal erosion and the enhanced BSA has subsequently been proposed.

Afterwards, the proposed procedure was applied to identify the parameters from the erosion coupled with mechanical tests of Hong Kong completely decomposed granite mixture (HK-CDG).

A coupled hydro-mechanical model of internal erosion

Brief introduction of the model

To assess the hydro-mechanical responses of a soil subjected to internal erosion, a novel hydro-mechanical model for internal erosion in granular soils developed by Yang et al.

(2019c) was adopted (summarised in Appendix). The saturated porous medium is considered to be composed of a stable fabric of the solid skeleton, erodible fines, fluidized particles, and a pure fluid phase. The internal erosion is considered as the phase transition of fine particles from a solid-like phase (described as erodible fines) to a fluid-like phase (described as the fluidised particles). The hydraulic process of suffusion is described by the mass balance of the mixture system and the phase transition of the fine particles is evaluated by an erosion law. The primary unknowns are the soil skeleton displacement (   c x t ). The variable x and t are space and time variables, respectively

A nonlinear incremental constitutive model (Yin et al., 2018a) was used to compute the effective stresses. The model was enhanced by introducing a non-linear elasticity, a nonlinear stress dilatancy law, the critical state concept, and the interlocking effect. It was then extended by defining the position of the critical state line in the log ' ep  plane as a function of the fines content in order to unify the mechanical behaviour of a sand-silt mixture from silt to sand or sand to silt.

The hydro-mechanical coupling of internal erosion has two main features: (1) the loss of fine particles affects the interlocking by changing the void ratio e and the critical void ratio c e simultaneously. When a part of the fine particles is detached from the soil skeleton, the void ratio e increases and the critical void ratio c e evolves with the fines content c f , creating an imbalance of internal and external forces and, consequently, a deformation of the granular assembly. Besides, the increase of the void ratio can reduce the elastic moduli G and K , leading to an increase of the deformation of the granular assembly; (2) the coupling from mechanics to hydraulics is considered implicitly by introducing the volume deformation in the mass balance equations. The change of the void ratio either increases or decreases, induced by the volume deformation of the soil skeleton. It may also affect the permeability as well as the flow rate and, inversely, the erosion rate. Moreover, the erosion model has been formulated to be strongly dependent on the shear stress ratio p M  according to the experimental results [START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF], while the impact of the confining pressure on the erodibility of the fine particles remains unclear or unsolved. Note that since the results in the literature present many inconsistencies [START_REF] Bendahmane | Experimental parametric study of suffusion and backward erosion[END_REF][START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF]Ke and Takahashi 2014a;[START_REF] Luo | Hydro-mechanical experiments on suffusion under long-term large hydraulic heads[END_REF][START_REF] Moffat | A hydromechanical relation governing internal stability of cohesionless soil[END_REF], the impact of the confining pressure was neglected as a first approximation.

Summary of model parameters

The proposed model (summarised in Appendix) requires 22 material parameters as well as the initial state parameters ( 0 e , 0 k and 0 c f ); 0 e is the initial void ratio, 0 k is the initial hydraulic conductivity, and 0 c f is the initial fines content. Based on their physical meanings, the parameters required by this model can be divided into five groups, as shown in Table 1. The procedure for determining the initial values of the state variables and some of the soil parameter values is relatively straightforward. The "standard" way to determine the model parameters is described as follows:

(1) The elastic parameters ( 0 K , 0 G and n ) can be determined directly from an isotropic compression test by plotting the isotropic compression line in the log ' ep  plane.

(2) The critical state-related parameters (  ,  and c  ) can be determined from a triaxial test. ), although their physical meaning is clear, their determination needs a procedure of trial-error by curve fitting, since they have strong coupled effects on simulating the evolution of the eroded mass and the hydraulic conductivity during and after erosion, as well as on the mechanical behavior before, during and after erosion.

The friction angle

Overall, some of the parameters can be determined directly based on the measurements, e.g., the elastic parameters and the critical state-related parameters, while some other variables cannot be measured directly, e.g., the parameters related to the interlocking-effect and the erosion process. Therefore, parameter identification using an optimization method is proposed for a relevant calibration of the coupled problem.

Identification methodology

In general, the procedure of optimization consists of three parts: a) defining a norm in a finite or infinite dimensional space, (b) expressing the mathematical formulation of the error by measuring the difference between the numerical and experimental results (in case of a finite space dimension), and (c) minimizing the error. An optimization strategy for searching the minimum of the error function has been used in order to reach this minimum.

Error function

To make the error independent of the type of test and the number of measurement points, an advanced error function has been adopted [START_REF] Levasseur | Soil parameter identification using a genetic algorithm[END_REF]. The average difference between the measured and the simulated results is expressed in the normalized form of the least square formula,   where

x is a vector of the unknown model input parameters; N is the number of the measurement points. exp i U is the value corresponding to the ith measurement point; num i U is the value corresponding to the ith point of the numerical simulation.

The scale effects on the fitness between the experimental and the simulated results can be eliminated by this normalized formula. Additionally, the objective error calculated by this function is a dimensionless variable. Thus, any difference in error can be avoided for different objectives with different variables.

Enhanced backtracking search algorithm

MBSA-LS

The MBSA-LS is a population-based evolutionary algorithm proposed by Jin and Yin [START_REF] Jin | Enhancement of backtracking search algorithm for identifying soil parameters[END_REF] based on the original BSA and differential evolution algorithm. It contains seven processes: initialization, selection-I, mutation, crossover, selection-II, local search and selection-III. A detailed description of the BSA can be found in [START_REF] Civicioglu | Backtracking search optimization algorithm for numerical optimization problems[END_REF]. The MATLAB code of the basic BSA can be found in BSA (2020). Here, the MBSA-LS is briefly described and the flowchart is shown in Figure 1. Step 2: Selection-I. At the beginning of each iteration, the historic population oldP is redefined through Eq.( 4), and then the order of the individuals in oldP is randomly changed by Eq.( 5). where a and b are two uniformly distributed random numbers between 0 and 1, and permuting function is a random shuffling function applied to a randomly selected previous generation as the historical population.

Steps 3-4: Modified mutation and crossover. The mutant individuals are generated by the historical and current individuals, as shown in Eq.( 6), which is different from the mutation of the original BSA. Then, the crossover operator is conducted by Eq. [START_REF] Calvello | Selecting parameters to optimize in model calibration by inverse analysis[END_REF]. A binary integervalued matrix (map) of size N×D guides the crossover directions of the BSA algorithm. The value of map is controlled by the mix rate parameter that is the only control parameter that should be determined in BSA, the details of the procedure can be found in [START_REF] Civicioglu | Backtracking search optimization algorithm for numerical optimization problems[END_REF]. Through the introduced mutation, some individuals are devoted to improving the population diversity by learning from the historical population, while other individuals are focused on enhancing the convergence speed via learning from the best individual in the current population. Consequently, the performance of the algorithm is enhanced with consideration of both the exploration and the exploitation abilities.

    old best     
Step 5: Selection-II. At this step, the population of the next generation new i P is generated according to a greedy selection mechanism. As shown in Eq.( 9), Vi is accepted if it provides a better function value than Pi, considering the minimum problem.

    , if
, otherwise where jrand=randint (1, D) is an integer randomly chosen from 1 to D and newly generated for each i, D is the dimension of the problem; the crossover probability CR∈  0,1 corresponds roughly to the average fraction of the vector components that are inherited from the mutation vector; CR=0.9 was used in this study.

i i i new i i V f V f P P P        .
To avoid a rapid loss of diversity, a competition between parents and children with an elitism strategy was adopted to perform the selection III. In the selection, 10% of individuals with the highest fitness are selected from the parents and children to survive in the next generation.

The remainders were chosen by tournament selection from the mating pool composed of parents and children apart from the 10% individuals.

Framework and pseudo code of the enhanced MBSA-LS

Based on the descriptions above, the pseudo-code MBSA-LS is summarized in Algorithm 1.

It can be observed that the proposed MBSA-LS retains the simple structure of the basic BSA apart from the implementation of the local search and the selection III.

Algorithm 1. Pseudo-code of MBSA-LS

Initialize population size and maximum number of iterations;

2. Initialize the current and historical populations using Eqs.( 2) and ( 3) for P and oldP;

3. Evaluate the errors of all individuals in the current population P; 4. iter=0;

5. For iter=1:Max_iter 6. Perform selection-I using Eqs.( 4) and ( 5) to form the historical population oldP;

7. Find the best individual from the current population;

8. Perform the new mutation using Eq.( 6) to generate the trial population M; 9. Perform the crossover using Eq.( 7) to generate the trial population V;  and  ) have strong combined effects on the evolution of the eroded mass and the hydraulic conductivity during erosion, as well as on the mechanical behavior before, during and after erosion, making it more difficult to be measured directly. Therefore, in this study, the elastic parameters, the critical-state-based parameters, the fines-content-related parameters and 1  have been measured directly, while for identifying the remaining erosion-related parameters and the interlocking-effect-related parameters, a mono-objective framework with six different criteria has been developed as follows:

            1~4 Êrror
                                                                 

 

Error a e   are, respectively, the error in the values of the deviatoric stress q and the void ratio e between simulations and experimental results for the triaxial shearing test. All the errors mentioned above were calculated according to Eq. (1).

Figure 2 shows the mono-objective identification procedure which aims to find values of the interlocking-effect-related parameters and the erosion-related parameters with the best attainable fit between model predictions and corresponding observations. In this study, the procedure was based on two different codes: the code for the integration of the coupled hydro-mechanical model of internal erosion to simulate the objective tests (Yang et al. 2019c) and the MATLAB code for the optimization process to find the minimum value of the error function. The erosion tests were carried out under initial isotropic and anisotropic stresses, with the major principal stress parallel to the seepage direction. After a saturation stage, the consolidation of the specimen was realized by gradually increasing the confining stress. For the tests under an anisotropic stress condition, the vertical stress was gradually increased up to the target value after the specimen had been isotropically consolidated to a given stress state. After the stress state being applied to the specimen, the pressurised water supply and the water collection system were connected to the triaxial system. The hydraulic gradient, i , was increased in stages to its final value. During the internal erosion testing, the applied hydraulic gradient, the vertical displacement, the outflow rate, the eroded soil mass and the specimen deformations were continuously measured. Conventional drained triaxial compression tests were then carried out on the specimens having been subjected to internal erosion.

A series of tests on the gap-graded HK-CDG mixture samples [START_REF] Chang | Internal erosion and overtopping erosion of earth dams and landslide dams[END_REF][START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF]) was selected for simulations, as summarized in Table 2. Note that the impact of the confining pressure on the erodibility of fine particles remains unclear and has not been considered in the present study. The erosion process with confining pressures of 100 and 200 kPa were not simulated. The physical properties of the soil mixtures are summarized in Table 3. could be easily obtained.

The erosion related parameter 1 0.92   was determined directly, based on the measurements of the final fines content after erosion.

All directly determined parameters are summarized in Table 4. They were fixed in the following parameter identification process.

Optimization and results

To identify the parameters related to the interlocking effect and the erosion process, it is necessary know the results of triaxial tests on the initial and eroded specimens as well as the erosion tests. Test ① is a triaxial test on the initial sample (see Table 2). It was first selected as an objective for the optimization process. Afterwards, to capture the stress-ratio-dependent features observed in the experimental tests, it is necessary to know the results of erosion tests under different stress ratios and the corresponding triaxial tests on the eroded samples in the optimization process. In this case, Test ④ and ⑧were selected as objectives for the optimization process.

According to Eq. ( 12), the objective error was defined as the average of

  Êrror mi  ,   Êrror ki  ,   Êrror a i   ,   Êrror r i  
for the erosion test and

  Error a q   ,

 

Error a e   for the triaxial shear tests. The initial population size was set at 160, and the number of generations was set at 80. The interval of each parameter is given in Table 5. Table 6 shows the optimal parameters and the corresponding objective errors. Using the optimal and the directly determined parameters, the objective tests were simulated with the adopted coupled hydro-mechanical model for internal erosion. Figure 7 shows the comparison of the eroded mass, the hydraulic conductivity, the axial and radial strains between the simulated and the experimental results of the erosion tests. Figure 8 shows the comparisons of the deviatoric stress and void ratio between the simulated and the experimental results of the triaxial tests. A good agreement between the simulations and experimental results was observed. Meanwhile, the parameters calibrated by a standard "trial and error" method (Yang et al. 2019c), and the corresponding objective error are listed in Table 6. It can be seen that the optimization method provides a set of optimal parameters with a smaller objective error, compared to the parameters calibrated by a standard "trial and error" method.

To evaluate the performance of the identified parameters, we used the optimal parameters and the directly determined parameters to predict the remaining erosion tests and the triaxial tests i.e., Test ②, ③, ⑤, ⑥, ⑦ and ⑨. Figure 9 and Figure 10 show the comparisons between predicted and experimental results of the erosion tests and the triaxial tests. In general, a good agreement between the simulations and observations indicates that the optimal values of the parameters are reasonable and reliable, which also suggests that the proposed identification procedure by the optimization method performs well in identifying the model parameters in the coupled hydro-mechanical model for internal erosion in granular soils.

As for some discrepancy between the predicted and the experimental results, this may be caused by deficiencies in the experimental tests and the inhomogeneity along the sample. In terms of the deficiency of the experimental tests, the lack of the time history of the hydraulic gradient for each test makes it difficult to identify the parameter related to the erosion rate.

Moreover, the number of experimental tests being limited, the selection of tests for identification of the parameters needs to be verified by conducting further experimental tests.

Another aspect which was not taken into account was the unknown influence of the inhomogeneity along the sample, which may lead to the reduction of hydraulic conductivity.

Not enough detail was given to calibrate the filtration term in the erosion law; the filtration was therefore not considered in the simulations. Furthermore, for the erosion test under a high stress ratio, especially when the stress ratio was close to the critical state, a small loss of fines content could change the soil from stable to unstable; the sample could become drastically inhomogeneous because of strain localization, which made it more complicated to predict the erosion-filtration process as well as the mechanical behavior after erosion.

Overall, using the optimization method can efficiently increase the accuracy of identifying the model parameters while minimizing the experimental costs and efforts, i.e., the "trialerror" efforts. Furthermore, the parameters determined by "trial-error" require some subjective judgment, and consequently this process is usually conditioned by the user's experience rather than being the pure product of some mathematical criteria. The use of optimization-based parameter identification can avoid such a problem. Therefore, the proposed identification of model parameters can be carried out by less-experienced operators, which should gradually promote the application of advanced hydro-mechanical models by the engineers.

Conclusions

This study has proposed an efficient optimization-based procedure for identifying the parameters related to soil internal erosion and interlocking effect using an enhanced backtracking search algorithm (MBSA-LS). In the proposed identification procedure, a newly developed coupled hydro-mechanical model was adopted to assess internal erosion with interlocking effect for granular soils. To identify the interlocking-effect-related and the erosion-related parameters, a mono-objective framework with six different criteria was developed, in which the elastic parameters, the critical-state-based parameters, the finescontent-related parameters, and the erosion related parameter were determined from experiments and kept constant in the optimization process.

From the erosion combined with mechanical tests of Hong Kong completely decomposed granite mixture (HK-CDG), the proposed procedure was then validated by identifying the parameters related to erosion and interlocking effects. The comparison between simulation and experiment results showed that the obtained optimal parameters were accurate and that their values were physically reasonable. The optimization method is shown to efficiently increase the accuracy of the identification of model parameters while minimizing experimental costs, as well as "trial-error" efforts.

Altogether, all results demonstrated that the identification of model parameters could be carried out by less-experienced operators, promoting the application of advanced numerical models in engineering practices. pressure. The model parameters are summarized in Table 1.

Figure captions

nR

  Constant controlling the degree of interlocking d Constant controlling the degree of interlocking Critical-state related  Constant controlling the slope of CSL  Constant controlling the nonlinearity of CSL Constant controlling the evolution rate of the transition zone d Ratio of the mean size of the coarse grains to that of the fine grains

cR

  is determined by plotting the stress path in the ' qp  plane,  and are determined in the log ' ep  plane according to the position of the critical state line.(3) the fines content related parameters, critical void ratios for the pure sand and the pure silt, which can be measured through critical states in is the ratio of the mean size of the coarse grains D50 to the mean size of the fine grains d50; a , erosion related parameter 1  controls the final residual fines content after erosion. It can be determined by comparing the final fines content to the initial fines content.For the remaining erosion-related parameters ( e

  .......................................................................... (1)

Figure 1

 1 Figure 1 Flowchart of MBSA-LS

P

  ................................................................................................. (9) Then, Steps 2-5 are continually performed until the termination criterion is satisfied. Step 6: DE-based local search. To further refine the quality of the current population, an elite method using the differential evolution as a local search has been introduced. The DEbased local search aims to find a better solution around the current best individual, which is performed after the BSA stage. The mutation strategy is presented as follows, .......................................................... (10) where the indices r1, r2, r3, and r4 are distinct integers uniformly chosen from the set is the best individual in the current generation i, which is randomly selected as one of the top 10% individuals in the current population obtained by Selection-II of BSA. The DE-based local search presented in Eq.(10) is employed to .............................................................. (11)

  10. Check the boundary for each trial individual 11. Evaluate the errors of trial population V; 12. Implement greedy selection-II using Eq.(9) to remain the better individual new P ; 13. Perform DE-based local search to generate new offspring DE new P ; 14. Perform the competition between DE new P and P to generate the new P for the next iteration; 15. End For3.3. Identification procedureThe elastic parameters ( , and ), the critical-state-based parameters( , and ) and 0 erosion related parameter 1  can be determined by comparing the final fines content to the initial fines content. However, the interlocking-effect-

    ......................... (12)where î is the hydraulic gradient; the errors of the eroded mass m  , the hydraulic conductivity k , the axial strain a  and the radial strain r  between simulations and experimental results for erosion tests;

FigureFigure 3 .

 3 Figure 2 Identification procedure flowchart



   was calibrated by plotting the stress paths of the drained triaxial tests in the ' pq  plane, as shown in Figure5. The critical-state-related parameters ( , ) were calibrated by plotting the critical state line in the log ' ep  plane measured from the drained triaxial tests before and after erosion (Tests ①-⑥). Along this way, curve fitting (see Figure4). The location of the critical state line in the plane is controlled by the initial critical void ratio 0 cr e which varies with the current fines content. The initial critical void ratio for the pure coarse and pure fine materials,

Figure 4 planeFigure 5 Figure 6

 456 Figure 4 Isotropic compression curve and critical state line in the log ' ep  plane

Figure 7 Figure 8 Figure 9 Figure 10

 78910 Figure 7 Comparisons between simulated and objective experimental results of erosion tests (a) Cumulative eroded soil mass; (b) hydraulic conductivity; (c) axial strain and (d) radial strain

  ................................................................................... are the porosity, the amount of erodible fines and the concentration of the fluidised particles, respectively; n is the source term describing the exchange between the erodible fines and the fluidised particles; s v is the velocity of the soil skeleton; ................................................................................................ (A5) with the permeability the medium k and the density of the mixture   ..................................................................................................... (A6) ........................................................................................................... (A7) where k  is the kinematic viscosity of the fluid; w p is the pore fluid pressure; s  and f  are the density of the solid and the fluid, respectively; k m is a power index which varies with the pore geometry.The exchange term n in Eqs. (A1)-(A4) is given as : .............................................................................................. (A8) ............................... (A9) where c f  is the ultimate fines content after a long seepage period; 0 c f is the initial fines content;  is the stress ratio; 1~3  ,  and e  are material parameters. The mechanical behaviour of the porous solid is controlled by: ............................................................................................................ (A10) where w p is the pore fluid pressure;  σ is the effective stress tensor, calculated by an fc- dependent non-linear incremental constitutive model (Yin et al. 2018a). The constitutive equations under axisymmetric triaxial condition are briefly summarised as follows: ........................................................ ................................................................................... .................................................. ............................................................................... ............................................................... (A15) with the critical state void ratio .......................................................................................................... (A16) ratio =' qp  is defined by the deviatoric stress q and the mean effective stress ' p ; d  and v are the increment of deviatoric strain and volumetric strain; pressure chosen as the reference stress; K and G are the bulk and shear moduli; 0 cr e is the reference critical void ratio at extremely low confining

Figure 1 Figure 2 Figure 4 planFigure 5 Figure 6 Figure 7 Figure 8

 1245678 Figure 1 Flowchart of MBSA-LS

Figure 9

 9 Figure 9 Comparisons between predicted and experimental results of erosion tests (a) Cumulative eroded soil mass; (b) hydraulic conductivity; (c) axial strain and (d) radial strain

Table 1

 1 Model parameters

  ................................................................................. (4)

	P	     	P	,	P	 , otherwise if 	, :	  0,1	
							
	old old	a b a b U permuting old old

 PP ....................................................................................................... (5)

  [START_REF] Nama | Improved backtracking search algorithm for pseudo dynamic active earth pressure on retaining wall supporting c-Ф backfill[END_REF] has been adopted to improve the performance of crossover as ................................................................................................ (8)

	follows,
	 mixrate 0.5 1   	rand	  0,1	
	M P F	 P P F p	
	     		, if , if	ij ij =1 =0
	a uniformly distributed random number within 	 0.3,1.0 controlling the amplitude of the
	search direction matrix, which is different from the F using in the basic BSA (Civicioglu
	2013) and the F in other similar algorithms that consider the best individual in the mutation
	ij	ij	ij
	V	P M		map map

P .............................................................................. (6) ............................................................................................................. (7)

where M is the population generated after mutation; pbest is the current best individual; F is

(Chen et al. 2017a; Chen et al. 2017b;[START_REF] Nama | Improved backtracking search algorithm for pseudo dynamic active earth pressure on retaining wall supporting c-Ф backfill[END_REF]

. Vij is the value of the jth variable for the ith trial individual. Furthermore, the mix rate of the ith individual proposed by

Nama et al. (

Table 2

 2 Summary of the testing program * ) CD1 and CD2 represent the drained triaxial test on the initial and eroded samples, respectively; ET represents the erosion test.

	Specimen identity	Test type ( * )	Confining Pressure (kPa)	Applied stress state Effective Mean Stress Deviatoric Stress (kPa) (kPa)	Stress ratio
	①BM-C-1	CD1	50	50	0	0
	②BM-C-2	CD1	100	100	0	0
	③BM-C-3	CD1	200	200	0	0
	④GS-C-1	ET+CD2	50	50	0	0
	⑤GS-C-2	CD2	100	100	0	0

(

Table 3

 3 Physical properties of the tested soil mixtures

Table 4 .

 4 Directly determined parameters of HK-CDG mixture

	 

Table 5

 5 Search domain and intervals of parameters related to interlocking-effect and erosion

		Interlocking-effect-related		Erosion related
			parameters			parameters
		Ad	np	nd	mk	e 2 3	
	Lower bound	0	1	1	1	1	1	1	0
	Upper bound	2	5	15	5 100 20 3	1
	Step Size	0.01	0.01	0.1	0.1 0.1 0.1 1 0.01

Table 6

 6 Optimal parameters with the errors * ) CD and ET represent the drained triaxial test and the erosion test, respectively.

	19 / 29

(

Appendix: Hydro-mechanical model of internal erosion (Yang et al. (2019)

  ...................................................................................................... (A1) ......................................................................................... (A2) .................................................................. (A3)

	A system of non-linear partial differential equations is deduced from the mass balance of the
	mixture system:
	  ws    div qv ...  0   div div  div        vv ...        div div         vv ...      div       
	cc	c s	s	c	s
	ss
	n cn f cc t ff tt tt	f	n

Acknowledgement

The financial supports provided by a GRF project (Grant No. 15209119) from Research Grants Council (RGC) of Hong Kong and the National Institute for Industrial Environment and Risks of France (INERIS) are gratefully acknowledged.