Numerical Study of Desiccation Cracking In Clayey Soil: A New FEM-MPM Coupling Method.

Desiccation cracking is a critical phenomenon soliciting the soil hydro-mechanical behavior, and significantly affects the performance of soil in geotechnical engineering. For this reason, an increasing interest toward studying and simulating the soil crack propagation, after a severe exposure to dry conditions (induced by desiccation), has been noticed during the recent years. However, major gaps remain in the previously developed models to properly provide a realistic prediction of the cracks pattern scheme especially when using the classical Finite Element Method (FEM), widely used in the geotechnical application. In this study, owing to the limitation of this method in re-meshing and dealing with large deformation, the authors were prompted to couple FEM with a mesh-free method: The Material Point Method (MPM) to overcome the individual drawbacks of each method. The dominant influencing factors on soil desiccation cracking have been assessed through a desiccation test performed in climatic chamber and using a digital image processing technique (image analysis) for a quantitative description of the studied sample. A model that relates porosity with suction and tensile strength was used to study the effect of the shrinkage phenomena in desiccation term, and to simulate the crack propagation in a thin clayey soil layer using the Code_Bright software. Consequently, a clear and connected crack pattern was observed, the problem of mesh dependency was clearly overcome proving the validity of the approach and providing a further insight into the behavior of clayey soil exposed to desiccation factors.


INTRODUCTION
The global warming effect has long been a question of great interest since it encompasses not only rising average temperature but also extreme weather: rainfall and droughts of unusual duration. As a consequence, soils undergo excessive swelling and shrinkage cycles, which can be a driving force for desiccation cracking [1,2].
The occurrence of soil desiccation cracking has a significant effect on structure safety as it causes crucial implication in the performance of soil such as leakage of nuclides in compacted bentonite [3][4][5][6], landslides and mudslides in clayey soil [7,8]. Understanding the phenomena is also of significant importance in the field of eco-geo construction such as raw earth building [9] which, despite being a prominent new technique, it faces the risk of desiccation cracking because of severe exposure to dry conditions, [10] especially during preparation of the material. Therefore, it becomes more important to study the behavior of clayey soil, and in particular the propagation of tensile cracks in earth structures.
As far as the mechanism of soil desiccation cracking is considered, cracks in Mode I are likely to occur when soil tends to shrink due to water loss from the soil surface to the atmosphere due to severe exposure to dry conditions. As soil in nature has a heterogeneous texture and a complex boundary conditions (compaction state, thickness, layer by layer decomposition, particle friction) shrinkage is restricted, internal tensile stresses develop and cracks initiate once the tensile stresses exceed the tensile strength of the soil [11][12][13].
Since the mid-twentieth century, multiple methods have been addressed to study and solve the issue. Typically, most of the published works concerned the laboratory and in situ observation of cracks formed by desiccation [3,[14][15][16]. Due to the limitation of the experimental approaches, various authors added the image processing technique to the experimental work for a better characterization of cracking pattern [10,[17][18][19] for a quick measurement of soil shrinkage characteristic curve. However, these authors focused mainly on factors influencing the process (drying conditions, specimen properties, boundary conditions etc.) rather than modelling and predicting. Therefore, other researches have made an effort by establishing numerical and theoretical schemes to simulate the evolution of crack pattern and clarify the mechanism behind cracks evolution during drying time. Finite element method (FEM) is one of the numerical methods that have been excessively used to simulate desiccation crack evolution [20][21][22][23]. However, being based on the continuum mechanics approach, FEM needed special treatment to properly simulate the cracking process noting the cohesive cracking mechanism and interface element method [24][25][26].
The finite element method is widely used in the treatment and simulation of fracture mechanic problems. However, the FEM can diverge when the size of one of the elements is equal to zero. Indeed, the interpolation functions are built on the mesh (Lagrangian, Eulerian). For this reason, one finds problems of mesh distortion. This problem is also present during the simulation of simultaneous evolving cracks.
To overcome the shortcoming of mesh based method like FEM, finite difference (FDM), finite volume method (FVM) etc., alternative spatial discretization methods broadly categorized into "mesh-free" and ''Particlebased" methods have also been developed to simulate desiccation induced cracking in soil [27]. For instance, [24,28] used the discrete element method (DEM) to represent the interaction between colliding particle during the shrinkage cracking process. Each component of soil is presented by a particle. The method shows a success in granular fields especially in clayey soils, though it is generally suitable for reduced-scale problems due to its high computational cost. Meshless smoothed -particle Hydrodynamics (SPH) was used by [29,30]. To model desiccation cracking, continuum equations are solved through spherical material points representing grains and fluid. (SPH) does not require any background mesh. Also, it is considered as the best suited method for large deformation and geomaterials failure problems.
The Material point method (MPM) of [31] considered as a non-dependent mesh method was used for the saturated and unsaturated deformable porous media in the geotechnical field. In fact, MPM combine the advantages of both Lagrangian and Eulerian description [27]. The media is discretized into multiple Lagrangian particles moving through a Eulerian mesh. Furthermore, unlike DEM, MPM allows the use of fully coupled hydromechanical process into the multiphase problems. Even though considered as a promising approach MPM usually have less computational efficiency than mesh based method, since Mesh-free interpolation and numerical integration require more computational resources. For these particular reasons, a number of algorithms that combine both approaches have been proposed [32] for meshless and FEM coupling technique. Besides, RPIM-FEM coupled algorithm showed a great potential in the work of [33,34].
This work aims at developing a continuum-based particle Method to capture the development of desiccation crack and overcome the shortcomings of existing numerical approaches. The developed method inherits all advantages of continuum mesh-based method of FEM and the strong point of the particle based method MPM. To fully understand the dynamic of crack propagation a desiccation experimental test was also carried out to observe the evolution of cracking pattern along with digital image processing technique for a quantitative description. The obtained experimental data is used to calibrate the model.
A Hydro-mechanical constitutive model by [35] was implemented in the finite element code CODE_BRIGHT software [36,37] to achieve the numerical analysis.

2-1. Soil characteristics
The Tibar-Beja clay, extracted from an area that has already experienced a landslide in November 2011, was used for this investigation. This type of soil has been extensively studied in the laboratory [38]. The geotechnical properties of the clay along with mineralogical composition and chemical analyses are summarized in Table 1.
The clay fraction is dominated by smectite and illite minerals. According to the unified soil classification system, it is classified as a high-plastic inorganic clay (CH).

2-2. Specimen preparation
The extracted soil was open-air dried and sieved through 2mm sieve. Over-saturated slurry was prepared by hand, mixing the Tibar-clay with distilled water in a plastic container at a gravimetric water content of about 130 %. After placing the plastic box in vibration device, the air bubbles are removed and clay mixture is left to dry for 3 days. A thin layer of 0.5 cm of water appeared above the slurry surface. Once pumped away a water content of 53% was measured.
For the desiccation test procedure the slurry was poured into rectangular (15 × 15 × 5 ) Plexiglas container with a 1cm layer thickness and then dried in a climatic chamber under a 35° degree temperature for 2 days in a row.

2-3. Test method
The drying of the prepared specimen was conducted in designed climatic chamber with imposed conditions (temperature T=35°c and relative humidity RH=40%). The chamber included a digital camera to automatically capture images at prescribed time frequencies and a sensor to capture the temperature and humidity inside the chamber.
The dynamic process of the crack propagation was recorded by the digital camera. Specimen was weighed every 2 hours to monitor the changes of water content during drying. Desiccation test was stopped when no change in specimen's mass was detected.

2-4. Image analysis:
A series of different crack patterns at different water contents was obtained. Results are illustrated below in Fig.  2. For a better exploitation of experimental data, extensive use is made of an image processing technique in order to quantify geometrical parameters of crack pattern. This method has been lately a success in the characterization of soil properties as a non-destructive method. [39] integrated X-ray computed tomography (CT) and digital image processing technique to investigate the evolution of desiccation crack network, [40] estimated the pore-size distribution of soil by applying tomographic and synthetic 3D images. [41] studied the effect of boundary conditions and environment conditions on soils cracking mechanism using image analysis technique.
In this work the specialized Crack Image Processing System (CIAS) was used to quantitatively describe the crack development. A set of geometric parameters was obtained following the image processing technique in (Table  2). The acquisition procedure is illustrated in six main steps: binarization (original image is transformed into black and white), noise elimination (black spots due to errors in quantification are removed), skeletonization (middle line of the crack segment is extracted from crack network to determine crack nodes and length), crack identification, and finally clod identification (see Fig. 1). More details about the image processing procedure using (CIAS) are presented in the manuscript of [18]. Fig. 1 Steps of image processing using (CIAS).
The topology of the images was used to quantitatively describe the characteristics of soil surface cracking [42]. The essential data to be obtained were crack ratio, total length of cracks, number of nodes along with the number of crack segments and number of clods in soil under desiccation conditions. The following measurements were carried out through the experiment test: gravimetric water content, suction, the state of saturation during cracking as well as the evolution of these variables and their possible correlation to cracking.   . 2 illustrates the crack initiation and propagation during desiccation process. It was found that at water content of about 25%, the first crack appears at the corner of the container, forming a 45° angle with the side walls which makes it in a good agreement with the work of [23]. Reporting that the first cracks are usually boundary cracks that initiate at the interface between the soil mass and the container wall. With further decrease in water content, we notice the propagation of "primary cracks" separating the surface slurry into several clods. With elapsed drying time, some new cracks known as "secondary cracks" initiate in clods or from already existing primary cracks. Fig. 2(d) indicates that at w=16.97% cracks intersection form a "T" shape and clods are divided into smaller ones, which is in a good agreement with the experimental results of; [5,18] Furthermore, it can be seen that at a very low water content w=5% already existing cracks gained larger width. It should be noted that crack width can continue to increase for a long period even after crack evolution come to an end [19]. In addition, at that final stage of drying, new cracks characterized by smaller width and length appeared from already existing primary and secondary cracks.

Fig
The Geometrical parameters obtained from slice images ( Table 2) provided a new insight into the cracking network within the soil. One can cite, the Crack intensity factor (CIF) defined as the ratio of the total crack area to total surface area [43]. In fact, Comparing the Crack intensity factor (CIF) evolution with decreasing water content in figure 4 highlights the value of water content corresponding to the appearance of the first crack w=25.6%. Accordingly, the evolution of CIF with elapsed drying time corresponds to the rapid decrease of water content. This is explained by the appearance of new cracks creating new boundary conditions that are in contact with the environment and which generate changes in the hydraulic boundary conditions in terms of suction. As a response to the matric suction increase, soil water content decreases [2]. As the water content approach to the shrinkage limit of the soil (w=15%), the crack development tends to slow. It should be noted that as drying time elapses, crack propagation gradually stabilizes and finally stops at CIF=16.7% corresponding to a water content of w=5%. The rest of geometric parameters (Crack area, total crack length, average crack width) follow the same trend as CIF while varying with decreasing water content. See for instance Fig. 3. To simulate the soil cracks induced by drying shrinkage, a hydro-mechanical model (in Section 3) was implemented in the Code_Bright software. Accordingly, in order to maintain the accuracy of the proposed approach all the hydraulic constitutive models were calibrated against the experimental test mentioned above. In fact, as the model is formulated in a net stress-suction framework, the soil water retention curve(SWRC) as well as suction variation with time are distinguished as the primary controlling factors of desiccation cracking [5]. Herein, we propose to calculate the soil degree of saturation using the geometry of resulting crack patterns experiments CIF and the soil mass variation during drying time as follows: Where is water mass; is the sample total surface; is the variating specimen height ; is the total solid volume. We should mention that changes of the layer thickness were monitored by a fixed camera as well during the experimental test. A linear decrease upon drying time was clearly observed. By further applying the Van-Genuchten Equation [44] suction variation with drying time can be obtained following the relation: where 0 is the capillary pressure parameter; λ is the shape's parameter and is defined as follows: where and are respectively the residual and the maximum saturation. = 0.05 ; = 1.
The driving idea behind the measurement technique of saturation degree using CIF and the image processing technique was to minimize the standard experimental measurement errors that neglect the crack evolution and effect during drying on soil properties. The Calibrated Suction increment during drying time (Figure. 5) was used as input data in the numerical model. According to [19,45] the soil suction is one of the main factors that influence the development of tensile stress and crack initiation and even used as a state variable for analysis [46]. For these reasons and for an accurate predicting of crack propagation scheme experimental calculated suction was used in the model.    Upon drying, the degree of saturation is related to the suction through the water retention curve (WRC) fitted in this work to the Van-Genuchten model (Fig. 6). In general, the degree of saturation of the soil decreases as matric suction increases. It is observed that the degree of saturation at low suction ranging from (0.01 to 0.1 MPa) the soil kept its stability (S r =1). Multiple authors have cited in their works that desiccation cracks initiate at soil suction close to the air entry value, which makes a good agreement with the experimental values in the water retention curve (WRC) (Fig. 6). Fig. 7 presents the decrease of saturation degree of the specimen with elapsed drying time. It can be observed that there is a slight difference betweem experimental result and numerical one especially from t=0.17 day, but generally, the curves follow the same trend. The numerical model seems to predict a rapid decrease of the saturation degree with a shorter drying time than the experimental one. It can be interferred that the numerical model is the first one to reach stablization before t=1day.

Numerical simulation
In this work, Simulation was performed to reproduce the dynamic development of the crack networkin soil and the trend of shrinkage due to desiccation and drying in general. According to [4,11] desiccation cracking is triggered by development of substantial matric suction in pore structure of granular soil due to drying phenomena. As a reaction, the soil will develop full potential shrinkage inducing an increment of horizontal tensile stress, followed by generation of vertical cracks when a particular tensile strength is reached. In fact, this process is very complex in nature because of the interactions that take place: the desiccation mechanism involves the interaction between both deformation of soil skeleton and development of pore-water pressure. For this particular reason and to introduce the non-linearity present during desiccation, a hydro-mechanical coupling approach is developed and implemented in the finite element CODE_BRIGHT software [36,37]. Thus, the implemented model is based on a set of Governing equations and constitutive laws replicating the correlation between the different hydraulic and mechanical properties (i.e. matric suction, tensile strength, cohesion, porosity and water content) that trigger the desiccation cracking development and propagation in general.

 Theoretical background and governing equation
The abundance of numerical approaches undertaken on desiccation cracking so far can be categorized into tensile stress failure approach and fracture mechanics approach [47]. Accordingly, a fracture phenomenon can either be depicted by a continuum formulation or a discrete notch.
The discrete approach found in multiple cases of work [26] simulate crack growth propagation by implementing cohesive joint between finite elements (dry bulks elements boundary). This denotes that the crack path has to be known in advance and undermines the ability to predict arbitrary crack propagation. Otherwise, if cohesive elements were implemented in all the internal elements' boundary an excessive increase in the degrees of freedom of the studied structure as well as the global effective stiffness will take place. Therefore, owing to their ability to successfully handle key issues associated to fracture development involving mode I fracture, the continuum failure criterion was suggested in this work.
As the soil is regarded to contain a set of interconnected pores that are randomly distributed, we assume that cracking is based on the simulation of the porosity rather than real discontinuity of cracks or fractures; this allows the simulation to maintain a continuous model rather than having to deal with discontinuity, which obviously offers certain advantages and meanwhile distinguishes itself from conventional fracture mechanics approaches such as the Extended Finite Element Method (XFEM), Generalized Finite Element Method (GFEM). Furthermore, Porosity appears in most of the governing equations (balance equation) conducting the multiphase flow of water and gas through porous deformable media like clayey soil. starting from the Mass balance of solid: Where ∅ is the porosity, and are respectively the mass and flux of solid. The variation of porosity caused by volumetric deformation is expressed through Equation (4) as follows: Where is the displacement. Porosity also appears in the equation of water mass balance Equation (6) not only as a function but also as a variation triggered by different processes, which clearly shows the influence of porosity variation in the balance equation of water.
Porosity plays an important role in the intrinsic permeability variation as well as the saturation degree of soil which clearly shows the influence of porosity in the balance equations and in the coupled flow-deformation process in general. Additionally [48,49] assumed that porosity has a significant effect on tensile strength and eventually, on cohesion of granular particles, in agreement with the experimental results of [35] as the direct tension test on specimen with different dry densities shows a dramatical decrease of tensile strength as porosity increases. This is explained by the broken links between soil particles .Moreover, in the study of [50] the cracks were considered as pore series of known geometry in the soils at different states. Eventually, porosity variation is considered as the keystone to successfully capture the development of fracture in clayey soil and used as an internal variable for the model. For the rest of the paper, our findings suggest that the cracked soil specimen is considered to be formed of two pore systems: cracks and soil pores. the two pore systems are characterized by different pore-size distribution. The pore size distribution for the soil matrix refers to the distribution of the total pore-volume in the soil matrix.

 Hydro mechanical model
To reproduce the main phases of desiccation cracking, a Hydro-mechanical framework is required. The coupling is performed through a sequential resolution of the two problems and the interaction between them. According to the literature reviews [51] desiccation cracking is developed essentially in Mode I fracture which clearly highlight the fact that fractures result from the progress of tensile stress between particles. As soil undergoes drying, it exhibits a suction change which entails a decrease of the water content and eventually an inter-particle force to increase. Consequently, an apparent cohesion term is thus created stimulating the emergence of tensile stress between dry bulk element.
The influence of suction and water content on soil mechanical characteristics (cohesion and friction angle) was already investigated over laboratory shear test [48,49]. The results emphasize the fact that the increase of suction in soil entails a decrease in moisture content and consequently an increase in cohesion and tensile strength.
To report the behavior of the desiccation process in a soil dominated by heterogeneity (initial distribution of porosity), the proposed model of [35] suggests the increase of the apparent cohesion under drying condition (i.e. suction rise) and the decrease of the same factor under porosity increase (Equation (7)). These assumptions are based on experimental tests established on the Tibar-Beja clay specimen (direct tension tests with different dry densities) and a proposed Modified Mohr coulomb criterion. For further details see the work of [35].
where is the cohesion and , are respectively suction variation and specimen porosity. The above formulation is introduced into mechanical formulation through elastoviscoplasticity. The formulation is based on the work of [52][53][54]. It allows the treatment of no-associated plasticity and strain softening behavior. Multiple authors [47,55] persisted on the importance of considering the elastoplastic behavior of soil during fracture.
The total strain rate ̇= is then given by the following relation: where the elastic behavior is carried out through a linear elastic formulation modelling the variation of young's modulus with suction term and porosity variation.
To reach the elastoplastic convergence, a regularization technique using the visco-plasticity was used. the idea is to soften the strain behavior by considering a degradation of the strength parameters when applying a viscosity term. This technique is only feasible by applying the following visco-plastic flow rule [53] : With is the fluidity parameter and denotes the relative rate of viscoplastic strain , is the stress function (it is also known as the flow function and increases monotonically with ) , is the yield function with the associative plasticity and G is the plastic potential. the stress function adopted is : correspond to the effective stress located outside the static yield surface ( > ) and is the power of the stress function which is an integer value and typically taken as 3. The inverse of viscosity can be written as a function of temperature : Where is considered as the activation energy. For temperature independent model, it is considered 0. The coefficient of viscosity(inverse of fluidity) is a crucial element in the viscoplasticity formulation: as the parameter takes a high values → ∞ the strain state will keep steady . This technique is known as the regularization technique. It is important to note that the viscoplastic strain occurs when stress state of the soil reaches a failure condition (where the static yield surface ≥ 0 ) The yield function and plastic potential are defined within the Drucker-Prager -based on Mohr-Coulomb failure criterion. For instance, to avoid the angular form of the coulomb yield surface the circular cone of Drucker and Prager has been chosen, which is perfect for low friction angle as in the case of our clayey soil.
In the plane (p, q') The yield function and plastic potential are defined by the following expressions: where p and q' are respectively the deviatoric stress and the effective mean stress. The above formulation allows a creation of shrinkage fracture in a consistent manner as if when the suction variation is low, the volume deformation is reversible and linear. However, when suction increases, the volume deformation shrinkage becomes significant which induces tensile stresses. Whenever the tensile stress reaches the tensile strength, porosity increases to reach reference porosity considered as the threshold for the crack emergence.
Moreover, it has been assumed that during the drying process, the fracturing behavior of the material is governed by net stress (total stress minus air pressure) and suction [56]. The effective stress could be then expressed in the following form [57]: where sigma denotes the total stress while and denotes respectively the air pressure and the water pore pressure while denotes the kroneker's delta It should be noted that this study focuses mainly on the desiccation crack initiation near the top surface of the soil. Thus, in view of simplicity, gravity forces are neglected in the governing equation of the problem. The model is implemented in Code-Bright software generating; a two dimensional flat layer of clay where the development of cracks is due to slow contraction of the material as a result of water evaporation [42].

 Evaluation of the model:
The ability of the proposed model to predict the soil desiccation cracking is further investigated in this section by simulating a restrained shrinkage induced soil cracking experiment. A square sample of clay set up initially as explained in the experimental test section is implemented in the Finite Element Code-Bright software. The Finite element mesh is discretized into 79 × 79 element representing the 15 cm × 15 cm soil slice area. Boundary conditions was induced in the computational mesh to impose specific constraints. As soil naturally possesses a heterogeneous fabric conditioned by mineralogical composition, a random field was generated numerically for a realistic simulation purpose. In fact, a Randomly distributed porosity is initiated with a given average initial porosity and a deviation by a 10 % of the given average value. The initial porosity is considered equal to 0.59. Porosity is then considered a function of time and space but with a constant values over each finite element. In the current simulation, the experimental suction profile is used as the shrinkage driver from which the effective stress is calculated. Fig. 8 shows the simulation results, using Code-Bright software, of the crack propagation presented in terms of porosity. When the tensile stress exceeds the tensile strength, cracks start to develop in the specimen. A complex network cracking can be observed representing the final state of drying obtained at t=1.025 days. The first four cracks appear simultaneously. These cracks start from the edge of the sample to the center dividing the specimen into smaller sub-soil cells. Cracks are seen to develop in both transversal and longitudinal direction. This can be explained by the inhomogeneous nature of the simulated specimen. In overall, the distribution of the cracks is not uniform across the specimen and crack pattern is less connected than a realistic network. To further verify the predictive capability of the proposed approach, model performance was evaluated graphically by comparing the experimental images obtained from the desiccation cracking laboratory test to the simulated images. Both converted to binary image, noise was removed, and crack edges were extracted. The area of cracks is measured by the total number of black pixels in the binary images. Image J program was used for the purpose. Simulated and measured Minkowski densities (CIF) versus time is presented in Fig .9. The comparison between simulation and experiments depict a very satisfactory agreement between both curve except for the slight difference in the interval between t= 0.7 and t= 0.89 days, and that could be attributed to the presence of imperfection in the experiment. Overall, the results show the capability of the model to capture and predict the response of the specimen to the desiccation conditions and the progressive development of cracks. The model seems to predict a shorter drying time than the experiment. After t=1.025 days crack evolution come to an end while in the experiment cracks propagation reaches stability at t=1.4 days.

FEM-MPM coupled method
Although the model provides a good description of the desiccation cracking process of clayey soil, a limitation in the obtained crack patterns was observed. Some cracks remain unconnected to each other and that is not frequent. This is probably because the approach is based on the finite elements. Every element behaves in an independent way. The (FEM) as a dominant numerical tool for fracture mechanics problems, often experiences such difficulties Furthermore, in mesh distortion and in simulating the propagation of simultaneous evolving cracks. For these particular reasons, alternative spatial discretization methods such as ''Mesh-less" and ''Particlebased" methods have been recently developed to overcome the drawbacks of mesh-based method such as (FEM, XFM (extended finite element) ….). For instance, Mesh-less and particle based methods have been receiving considerable attention lately, since this class of methods shows great potential to meet the demands of modern software, error estimators, though they usually have worse computational efficiency than (FEM)., since Mesh-free interpolation and numerical integration requires more computational resources.
In this paper, our findings suggest to develop a continuum-based particle Method to capture the development of desiccation crack and overcome the shortcoming of existing numerical approach. The Material Point Method (MPM) is one of the most distinguished particle method, its capacity to deal with large displacement and its natural dynamic formulation make this technique an important numerical tool to tackle a number of geotechnical problems [58]. For instance, the developed method will inherit all the advantages of a continuum mesh-based method of FEM and the strong point of the particle based method (MPM).
This study focuses mainly on the development of a (FEM-MPM) coupled method and its application to simulate the desiccation cracking in clayey soils and capture the process of crack growth.

 Material point method
The (MPM) of Sulusky [31] can be viewed as an extension of the classical (FEM) . Soil and structural bodies are represented by Lagrangian particles called Material points or particles that move through a Eulerian fixed mesh. The physical properties reside with particles throughout the computations ,whereas the Eulerian mesh carry no permanent information [59].
In fact, the first two stages of the two numerical methods (FEM and MPM) are similar. The similarities between the particles and the Gauss point are clearly observed when comparing figures in Fig. 10 although the difference reside in the last stage. Moreover, the Gauss point in the finite element method is fixed to the cell it was created in. However, in the material point method, particles are free to move from one cell to another to avoid the mesh interference. During the convective phase, the particles remain in their updated positions, while the mesh is reset to its initial one.
From the continuum point of view, the (MPM) particles represent discrete samples of the continuous material and the grid is just a helper for computing their physical interactions. In the work of [60][61][62], on a sandy soil specimen, the particles were rendered as individual grains. Herein, for the proposed numerical approach we propose the following assumption: porosity will be considered as the particle representation of the continuum. The driving idea behind this approach is in the work of [36] porosity is defined as an element-wise in the numerical scheme .In fact, porosity is considered a function of time and space but with constant values over each finite element.

 Proposed approach steps
The essential steps of the proposed approach continuum-based particle Method are summarized in the algorithm below: In the beginning of the analysis, the entire problem domain is discretized into finite elements for computational efficiency. Crack propagation is then generated by the CODE-BRIGHT software. Once the crack pattern is established and the step of the mesh divergence is detected, the approach of the (MPM) is then introduced: The continuum is discretized into porosity elements. The porosity elements are from now on defined as the new material points whose displacements are detected in each time increment. Once the mesh deforms in such a way that some elements become heavily distorted, leading to numerical inaccuracies and a potential breakdown of the calculation, the crack scheme right before that critical step is captured via an image analysis technique software and a created maple subprogram, elements are converted into pixels and then into porosity particles. In order to quantify crack pattern, the area of cracks is measured by the total number of black pixels in the binary images and the physical resolution(e.g. one pixel per 0.379 mm) of the image . It is defined as follows : each slice contains 395*395 pixels with a pixel size of a 0.379) . The created grid acts like a scratch pad and can be destroyed and reset after each solve, and reinitialized in the beginning of the next time step. The reconstructed crack scheme was then injected in the new computational framework and cracks could pursuit the propagation until another divergence is detected.
In literature, different numerical techniques have been used as enhancement for the (FEM) to deal with mesh distortion. These enhancements have managed to tackle these difficulties with different degrees of success. For example, the remeshing technique can be computationally expensive and also troublesome with history-dependent materials due to the mapping of material state variables from the old to the new mesh [26]. Also from that point of view, MPM avoids both the Lagrangian problem of mesh distortion and the Eulerian problem with the convective term: The Lagrangian material points or particles eliminate the convective term of the material time derivative, while the Eulerian background grid avoids mesh distortion. However, When particles move through the mesh, there is a jump in the acceleration of the particle which leads to a possible numerical inaccuracies in the mapping phase and unbalance in the nodal internal forces. For all this multiple reasons, we developed this new method in order to overcome both the drawbacks of remeshing FEM and MPM.

 Simulation of desiccation cracking using the proposed approach
The proposed approach is further applied in this section to simulate soil cracking problem due to desiccation phenomenon. Fig. 11 clearly displays the continuity of crack propagation through the proposed approach. In fact, the obtained crack scheme from the above FEM simulation (step1) suffering from crack discontinuity was treated at t=0.85 days when the simulation endures a mesh distortion. The crack scheme at the critical step is captured (step2) and was injected in a new one. The potential of the proposed simulation approach is demonstrated in (step 4) through the dichotomy process. In fact, a subsequent crack appears in the middle of two existing neighboring ones [2,11,63] and assure the continuity of developed cracks (step 4). The results in Fig. 11 show the ability of the new approach to simulate a realistic crack pattern. Crack propagation divided the specimen into new sub-cells which make it in a good agreement with the experimental observed crack pattern. At the end of the simulation (step 4), fully propagated cracks are formed which depict the capability of the developed approach to simulate and predict a crack pattern close to reality. Toward the end of the desiccation cracking simulation, the development of cracks tends to stabilize and no new cracks develop. It is important to note that when using the proposed approach, the type of element must also be chosen carefully in order to eliminate potential inaccuracies caused by discontinuities in element dimension. These consideration lead to the use of quadratic 4 nodes elements.
In general, element stresses and strains also need to be interpolated into the new generated mesh. However, since an elastoplastic material is assumed, the new values could be calculated directly from the input porosity. In addition, nodal masses are recalculated to correspond to the new mesh. This differentiates the method from the conventional remeshing strategies that require transferring a substantial amount of data from the old mesh to the new one. Results show this algorithm to be quite reliable. Fig.  12 depicts the considerable tensile stress decrease in crack tips as the particle porosity reaches the reference porosity This leads it to be in a good agreement with the proposed model of crack formation in mode I. Accordingly, tensile stresses locally develop, and when reaching the maximum strength, the soil fails. The crack is formed and tensile stress decreases. This is because energy is released at the crack tip. Fig. 12 Tensile stress evolution at a crack tip. Fig. 13 Plastic deformation evolution at a crack tip.
Strain variable are also generated in post processing. The results show in Fig. 13 Plastic Deformation generated at both cracked and non-cracked zone. The figure clearly depicts the large captured strain evolution around the crack tip.

CONCLUSION
Desiccation cracking in soil is common in geotechnical field and can cause many safety issues. An understanding and modelling of the cracking mechanism would enable mitigation and prevention of these issues. Even though the heterogeneity and nonlinearity of soils, as well as the fact that it is not a material made under factory conditions like many other engineering materials has been a factor in making modelling its exact behavior difficult. In this study, a continuum-based particle Method (FEM-MPM) coupling the advantages of both approaches, is developed and implemented to simulate the initiation and propagation of crack pattern in clayey soils subjected to desiccation. A Hydro-mechanical model including cohesion suction and porosity evolution is proposed as a constitutive equation. The simulation has been performed using the CODE_BRIGHT software.
To develop a more refined model of desiccation, cracking Laboratory test results were used as a basis for validation and calibration of the model. The proposed approach is verified by simulating a desiccation cracking test. The proposed framework managed to reproduce the fracture pattern observed in the experiments. Moreover, the proposed model features a clear and connected crack network: proof of the potential of the FEM-MPM coupled method, proposed in this paper to form nonlinear deformation under form of desiccation cracking. In future work the FEM-MPM is seen to simulate the progression of the landslide and generate a reasonable post-failure configuration of the problem. The present approach will be extended in further studies to include the effect of water in cracks. Figure 1 Steps of image processing using (CIAS).

Figure 2
Evolution of desiccation cracking on the surface of Tibar-Beja clay samples.

Figure 3
Variation of cracks geometric parameters with gravimetric water content.

Figure 4
Changes in water content and CIF with drying time.     Evolution of CIF (Experimental and numerical result).  Simulation of crack propagation between consecutive steps (Elements number 6241).

Figure 12
Tensile stress evolution at a crack tip.

Figure 13
Plastic deformation evolution at a crack tip.