A coupled hydro-mechanical model for subsurface erosion with analyses of soil piping and void formation

A coupled hydro-mechanical erosion model is presented that is used for studying soil piping and erosion void formation under practical, in-situ conditions. The continuum model treats the soil as a two-phase porous medium composed of a solid phase and a liquid phase, and accounts for its elasto-plastic deformation behaviour caused by frictional sliding and granular compaction. The kinetic law characterizing the erosion process is assumed to have a similar form as the type of threshold law typically used in interfacial erosion models. The numerical implementation of the coupled hydro-mechanical model is based on an incremental-iterative, staggered update scheme. A one-dimensional poro-elastic benchmark problem is used to study the basic features of the hydro-mechanical erosion model and validate its numerical implementation. This problem is further used to reveal the interplay between soil erosion and soil consolidation processes that occur under transient hydro-mechanical conditions, thereby identifying characteristic time scales of these processes for a sandy material. Subsequently, two practical case studies are considered that relate to a sewer system embedded in a sandy soil structure. The first case study treats soil piping caused by suffusion near a sewer system subjected to natural ground water flow, and the second case study considers the formation of a suffosion erosion void under strong ground water flow near a defect sewer pipe. The effects on the erosion profile and the soil deformation behaviour by plasticity phenomena are elucidated by comparing the computational results to those obtained by modelling the constitutive behaviour of the granular material as elastic. The results of this comparison study point out the importance of including an advanced elasto-plastic soil model in the numerical simulation of erosion-driven ground surface deformations and the consequent failure behaviour. The numerical analyses further illustrate that the model realistically predicts the size, location, and characteristic time scale of the generated soil piping and void erosion profiles. Hence, the modelling results may support the early detection of in-situ subsurface erosion phenomena from recorded ground surface deformations. Additionally, the computed erosion profiles may serve as input for a detailed analysis of the local, residual bearing capacity and stress redistribution of buried concrete pipe systems.


Introduction
Soil erosion by water is a degradation process in which particles are detached from the soil structure and carried away by water flow. This process can be activated both by overland and subsurface flow [3]. When soil erosion takes place by subsurface flow, the time and spatial development of the erosion profiles are not visible, so that the process may remain unnoticed for a long time. This aspect makes it generally difficult to monitor and analyse subsurface erosion, which is likely the reason that the number of research studies on subsurface erosion -also referred to as internal erosion -are disproportional compared to those on surface erosion [3,48]. However, subsurface erosion has been reported as a highly relevant and widespread process [48], and is the cause of various catastrophic failure phenomena with large societal impact, such as sewer collapse, dyke break-through, the formation of sinkholes and the instability of earth dams [13,19,24,67,69].
An important subsurface erosion process is soil piping, whereby relatively large porous networks are created under preferential flow paths [3,4,32]. Soil piping has been observed in various types of soils, such as loess-derived soils, organic soils, clayey soils, and sandy soils [4], and may be initiated when the soil contains fine particles that are smaller than the average pore space in the soil. When the process starts and small particles are eroded, the pore space increases and individual pores become connected and form channels, which enlarges the hydraulic permeability and can also trigger the release of coarser particles [3,32]. The pipe network that develops may lead to subsidence of the soil structure, and eventually even to severe structural failure, such as an abrupt collapse of the pipe network and the overlying soil structure -termed a sinkhole -, or landsliding [3,67]. Whether or not a pipe network collapses depends on the loading conditions and on the level of migration of soil particles. If, under the specific hydraulic conditions applied, only the finer particles erode and the coarser particles to some extent maintain a particle contact force fabric, the soil structure preserves load bearing capacity and does not collapse. This process is denoted as suffusion [22,30,43]. Conversely, the erosion process is characterized as suffosion if under a strong water flow the transport of fine particle is accompanied by the collapse of the contact force fabric of the coarser particles, whereby some form of catastrophic structural failure is induced [22,30,43].
Another important type of subsurface erosion is void formation, which may occur around buried concrete pipe systems, such as sewer pipes, industrial discharge lines, tunnels, culverts and storm drains [42,47,68]. Erosion voids can emerge due to ageing or improper construction procedures of pipe systems, with defects (i.e., cracks or gaps) in the pipe system initiating a water flow that locally washes away the soil support. Under continuous growth the erosion void may eventually evolve into a suffosion sinkhole, whereby the bearing resistance of the supporting soil vanishes and the soil structure above the pipe collapses [31,67]. Sinkhole failures have been reported to happen within the lifetime of the pipe system, as well as during the period of construction [67], indicating that the rate of the underlying erosion processes is strongly variable and determined by local hydro-mechanical conditions. Various experimental and numerical studies have addressed the effect of erosion voids on the stresses generated in pipe systems and in the surrounding soil, whereby the stress redistribution caused by the appearance of an erosion void is found to be characterized by its location, size and contact angle [41,47,61,68,71]. Since under void erosion the load transfer from the surrounding soil to the pipe generally takes place more localized, whereby the load magnitude increases, the susceptibility of the pipe system to cracking and damage, and thus to catastrophic failure, also increases [42,51,60,68].
The above examples illustrate the importance of accurately predicting subsurface erosion processes by physicsbased modelling. Erosion models can be generally divided into two categories, namely interfacial erosion models and bulk erosion models. In interfacial erosion models, the local erosion process at the soil-fluid interface of a water channel is considered to be a ''threshold phenomenon'', which initiates and develops once the critical soil erosion resistance is exceeded [29,34]. Since interfacial erosion is a shear induced mechanism, the rate at which the soil particles are released is typically described by the amount at which the flow shear stress at the soil-fluid interface exceeds a critical threshold value [1,12,11,16,34]. Interfacial erosion models have been applied to predict the amount of soil detachment and the characteristic erosion time for individual flow channels in different types of soil, and subjected to various flow conditions [10,11,34]. This has illustrated the complexity of comparing experimental erosion data and extending the results to conditions different from which they were acquired. In bulk erosion models -also referred to as continuum erosion models -the erosion of the complex, microstructural pore network of a soil is modelled as an effective, hydrological process at material point level, which satisfies the mass-balance equations and particle transport kinetics for a porous granular medium saturated with a fluid [5,17,50,63]. As an extension, the mechanical behaviour of the particle structure may be incorporated by coupling such a hydrological model with a constitutive soil model and accounting for the balance of linear momentum. In accordance with this principle, a limited number of coupled hydro-mechanical formulations have been proposed in the literature, which consider different constitutive formulations and various degrees of coupling [46,53,69]. By means of illustrative examples, these works clearly demonstrate the value of a coupled hydro-mechanical soil erosion model for the identification of the time instant at which erosion-driven structural failure starts, and the extent to which failure occurs. Hence, the erosion profiles computed by these models may be used as input for a detailed analysis of the local, residual bearing capacity and stress redistribution of buried concrete pipe systems [42,47,51,71]. Additionally, coupled hydro-mechanical erosion analyses may support the early detection of in-situ catastrophic soil erosion phenomena, as registered from temporal changes of soil surface profiles measured by satellite radar interferometry [15,40].
In order to explore these challenging aspects in detail, in the present communication a coupled hydro-mechanical bulk erosion model is presented that is used for studying soil piping and erosion void formation under practical, insitu conditions. The model treats the soil as a two-phase porous medium composed of a solid phase and a liquid phase, and accounts for its elasto-plastic deformation behaviour caused by frictional sliding and granular compaction. The kinetic law characterizing the erosion process is assumed to have a similar form as the type of threshold law used in interfacial erosion models. The degradation of the particle structure under erosion is considered to reduce the effective elastic stiffness of the porous medium and increase its permeability. The reduction in elastic stiffness is accounted for via an erosion degradation parameter, which depends on the actual porosity of the soil and quantifies how effectively the mechanical degradation under erosion takes place. The coupled hydro-mechanical model is implemented in the commercial finite element modelling program ABAQUS, whereby the solutions of specific problems are computed in an incremental-iterative fashion using a staggered update scheme. As a start, a onedimensional poro-elastic benchmark problem is considered, which demonstrates the basic features of the hydromechanical erosion model and validates its numerical implementation. This example is further used to reveal the interplay between soil erosion and soil consolidation processes that occur under transient hydro-mechanical conditions, thereby identifying characteristic time scales of these processes for a sandy material. Subsequently, two practical case studies are considered that relate to a sewer system embedded in a sandy soil structure. The first case study treats soil piping caused by suffusion near a sewer system subjected to natural ground water flow, and the second case study considers the formation of a suffosion erosion void under strong ground water flow near a defect sewer pipe. The effects on the erosion profile and the soil deformation behaviour by plasticity phenomena are elucidated by comparing the computational results to those whereby the constitutive behaviour of the granular material is modelled as elastic. The results of this comparison study illustrate the importance of including an advanced elasto-plastic soil model in the numerical simulation of erosion-driven ground surface deformations and the consequent failure behaviour.
The paper is organised as follows. In Sect. 2 the governing equations for the process of subsurface soil erosion are formulated. Section 3 demonstrates how the effect of erosion is incorporated in an advanced elasto-plastic model for the frictional sliding and granular compaction behaviour of the soil. The incremental formulation is based on the flow theory of plasticity and the assumption of small deformations. The incremental-iterative time integration scheme of the constitutive model is outlined, and the calibration of the material parameters is discussed for a sandtype material, using test results presented in the literature. In Sect. 4, the mass-balance equation for the pore fluid established in Sect. 2 is further developed into the governing differential equation for the pore pressure field, thereby identifying the modelling terms contributing to soil erosion and soil consolidation. Subsequently, the staggered solution procedure of the coupled hydro-mechanical erosion model is presented. By means of a one-dimensional benchmark problem, in Sect. 5 the basic features of the coupled erosion hydro-mechanical model are demonstrated and its numerical implementation is validated. In addition, the interaction between soil erosion and soil consolidation processes is considered through the identification of the corresponding characteristic time scales. Section 6 presents two practical case studies related to a concrete sewer pipe embedded in a sandy soil structure. The first case study treats the process of soil piping due to suffusion under natural groundwater flow, and the second case study considers void formation by suffosion under strong groundwater flow towards a defect sewer pipe, i.e., a sewer pipe with a gap created by a sudden failure of a pipe connection. Finally, Sect. 7 summarizes the main conclusions of the work and provides recommendations for future research.
As a general scheme of notation, scalars are written as lightface italic symbols and vectors and tensors are indicated as boldface symbols. The inner product between two vectors a and b is indicated by a centered dot as a Á b (i.e., a i b i ), and the action of a fourth-order tensor A on a secondorder tensor c by a double dot product A : c (i.e., A ijkl c kl ). Further, $ represents the gradient operator.

Subsurface erosion
This section formulates the governing equations for the process of subsurface soil erosion, which serve as input for the constitutive soil model presented in Sect. 3 and the coupled hydro-mechanical numerical formulation presented in Sect. 4. The erosion model is based on the assumption of a two-phase porous medium composed of a deformable soil skeleton saturated with pore fluid(s), for which the basic principles were originally established in the landmark contributions of Terzaghi [62] and Biot [6] for the study of soil consolidation. The porous media theory has been subsequently extended and generalized by many others for a range of applications, see the reference works [7,8,18,21] for an overview. The formulation below further relies on concepts and notions presented in [63] for the hydrological modelling of bulk erosion processes.

Governing equations
Consider a granular material of which the pore space is fully saturated with a fluidic mixture. The total elementary material volume dV and elementary mass dm of the twophase porous medium can be decomposed into a solid (s) part dV s with mass dm s and a fluidic mixture (fm) part dV fm with mass dm fm as In addition, the densities of the solid and fluidic mixture phases respectively follow from Dividing Eq. (1) 1 by the elementary material volume dV turns this expression into with the volume fractions of the solid and fluidic mixture phases given by / s ¼ dV s =dV and / fm ¼ dV fm =dV, respectively. Note from Eq. (3) that / fm may be alternatively referred to as the porosity of the porous medium. In accordance with the above definitions, the bulk densities per unit bulk volume for the solid and fluidic mixture phases are respectively given bỹ Further, with Eqs. (1) 2 and (4), the bulk density of the twophase porous medium reads During subsurface soil erosion, solid particles are released from the particle skeleton and carried away by the fluid flow in the pores, a process known as ''fluidization''. The transport of solid particles by the pore fluid causes the fluidic mixture introduced above to be composed of fluid (f) particles and fluidized solid (fs) particles, see also Fig. 1.
The velocities v f of the fluid particles and v fs of the fluidized solid particles may be assumed to be equal to each other, and thus equal to the effective velocity v fm of the fluidic mixture [63] When denoting the elementary mass and volume of the fluidized particles as dm fs and dV fs , and the elementary mass and volume of the fluid particles as dm f and dV f , the density q fm of the fluidic mixture can be developed from its definition in Eq. (2) as with c fs ¼ dV fs =dV fm the concentration of fluidized particles in the fluidic mixture and q f ¼ dm f =dV f the density of the fluid. For a sandy soil containing pore water, the values of q s and q f in Eq. (7) correspond to q s ¼ 2650 kg/m 3 (quartz) and q f ¼ 997 kg/m 3 (water). Note further that in the limits of c fs ¼ 0 and c fs ¼ 1 the density q fm of the fluidic mixture correctly reduces to the density q f of the fluid phase and the density q s of the solid phase, respectively.
In correspondence with the above definitions, for the process of soil erosion the mass-balance equations for the solid and fluidic mixture phases are respectively expressed by [38,63] in which v s is the velocity field associated to the solid phase and _ m er is the rate of net mass eroded and fluidized at an any time instant and location. Since _ m er [ 0, from Eq. (8) it is clear that the fluidic mixture phase receives mass from the solid phase. When the frame of reference used for describing the deformation of the two-phase porous medium is attached to the solid phase, the velocity v in a material point equals the velocity v s of the solid phase In addition, it is convenient to introduce a relative velocity between the fluidic mixture phase and solid phase as [38,63] v Inserting Eqs. (9) and (10) into Eq. (8) results in with the superimposed dot henceforth indicating the rate with respect to the frame of reference attached to the solid phase. The summation of Eqs. (11) 1 and (11) 2 leads to the total mass balance equation of the two-phase porous medium with the bulk densityq of the two-phase medium given by Eq. (5). Note that the left-hand side of Eq. (12) contains similar terms as the mass balance of the solid phase, see the left-hand side of Eq. (11) 1 , while the right-hand side represents the net mass discharge by the fluidic mixture present in the pores. Obviously, the eroded mass does not appear in the total mass balance, as, from the conservation of mass, there is no mass production/reduction in the assembly of the solid and fluidic mixture phases. Substituting Eq. (4) into Eq. (11) gives in which, for simplicity, both the solid and fluidic mixture phase materials are considered as incompressible, so that their densities remain constant in time, oq s =ot ¼ oq fm =ot ¼ 0. Note that the mass-balance equation of the solid phase, Eq. (13) 1 , determines the actual value of the solid phase fraction / s , which thus may serve as an internal state variable for the erosion process. The summation of Eqs. (13) 1 and (13) 2 leads to the inclusion of the effect of eroded mass in the so-called storage equation: in which Eq. (3) and its time derivative oð/ s þ / fm Þ=ot ¼ 0 have been applied, as well as the relation between the flux q fm of the fluidic mixture (also referred to as the Darcian velocity) and the relative velocity v r [38,52,63] For solving the mass-balance equations in the forms given by Eqs. (13) 1 and (14), the actual density q fm of the fluidic mixture present in the pores needs to be calculated. As indicated by Eq. (7), this density depends on the actual concentration c fs of fluidized particles in the fluidic mixture. In general, c fs can be computed by treating it as an internal state variable in a separate mass-balance law for the fluidized particles [63]. However, when the characteristic time scale of the erosion process is (much) larger than the characteristic time scale of the pore fluid flow process, the actual concentration of fluidized particles in the fluidic mixture typically is relatively small, c fs ( 1. Under this condition, which is applicable to the erosion problems analysed in the present work, see Sects. 5 and 6, it is reasonable to assume that the concentration c fs takes a (very) small, constant value. Correspondingly, the density q fm of the fluidic mixture remains constant, and follows from Eq. (7) by substituting the specific value of c fs . The rate of eroded mass _ m er appearing in Eqs. (13) 1 and (14) needs to be defined by means of a kinetic law. In the present work, this kinetic law is considered to have a similar form as the type of threshold law regularly used in erosion models at fluid-soil interfaces in channel and pipe flows, whereby the flow shear stress at the wall needs to exceed a critical value in order to generate mass transport of particles [1,12,16,34]. Since in bulk erosion models the interfaces between individual particles and the pore fluid are not modelled explicitly, the interfacial shear stress appearing in such threshold laws needs to be replaced by an appropriate, continuum-type of process parameter. For this purpose, consider the case of pressure-induced laminar flow of a Newtonian fluid through a pipe, i.e., Hagen-Poiseuille flow, whereby the average shear stress at the pipe wall is proportional to the average flow velocity in the pipe [37]. Hence, from the analogy between a pipe and a pore channel, the erosion kinetics of the granular medium may be considered to be governed by the relative fluid velocity v r given by Eq. (10). In accordance with the above argumentation and Eq. (15), the initiation and growth of erosion are formulated by expressing the rate of the eroded mass _ m er in terms of the flux q fm of the fluidic mixture as _ m er q s ¼ k er kq fm k À q fm;cr À Á / fm if kq fm k [ q fm;cr and / s ! / min s ; 0 otherwise : Here, k er is the erosion coefficient (with dimension length À1 ), kq fm k is the Euclidean norm of the flux of the fluidic mixture, q fm;cr is the critical flux threshold above which erosion takes place, and / min s is a prescribed, minimum value of the solid volume fraction below which erosion will not occur. Eq. (16) illustrates that the erosion process is active when i) the fluid flux exceeds the specific threshold value and ii) the solid volume fraction has not (yet) reached a prescribed, minimum value. Note from Eq. (16) that the erosion process is more intense when the porosity / fm is smaller, i.e., in the case of smaller pore canals. Further, the critical flux threshold q fm;cr is considered to be somewhat analogous to the so-called Shields number [44], which is a dimensionless, critical shear stress commonly used in interfacial erosion modelling for calculating the initiation of sediment motion in fluid flow. The critical value q fm;cr needs to be determined experimentally, and, as indicated by the Shields number, is expected to scale with the characteristic particle size of the sediment.
In accordance with Terzaghi's law [62], the total stress r experienced by the saturated granular medium may be decomposed into an effective stress r 0 in the particle skeleton and a pressure p fm in the fluidic mixture occupying the pores: where 1 is the second-order unit tensor. Following the sign convention typically used in the solid mechanics research community, the stress (pore pressure) is assumed to be positive for tension (suction) and negative for compression. Similarly, the volumetric deformation generated under the applied loading is assumed to be positive for dilation and negative for contraction. During the erosion process, the degradation of the granular skeleton transferring the effective stress r 0 in the present work is assumed to be characterized by an erosion degradation parameter D er of the form: In the absence of erosion, the solid volume fraction remains constant and thus is equal to its initial value, / s ¼ / s;0 , which sets the above degradation parameter to zero, D er ¼ 0. Conversely, a fully eroded material point relates to / s ¼ / min s ¼ 0, in correspondence with the degradation being equal to unity, D er ¼ D max er ¼ 1. Note that a maximal degradation parameter D max er lower than unity corresponds to a prescribed minimum solid volume fraction larger than zero, / min s [ 0; this condition is assumed to be representative of suffusion erosion, whereby the granular skeleton does not fully degrade and thus preserves load bearing capacity. Additionally, the exponent b determines how effectively the mechanical degradation of the granular skeleton takes place during erosion. Specifically, a well graded soil typically contains a substantial amount of fines that occupy the void space and do not contribute to the contact force fabric, so that the effect of erosion of these fines on the initial mechanical degradation of the particle contact force fabric is minor, in correspondence with a relatively high value of b. Conversely, for a poorly graded soil the particle contact force fabric is more easily degraded by the erosion process, in correspondence with a relatively low value of b. Due to a lack of experimental data, in the present study the parameter b is arbitrarily set equal to unity, b ¼ 1:0. In Sect. 3 the degradation parameter D er will be used to consistently include the effect of erosion in an elasto-plastic constitutive model for the soil.
The fluid flux q fm appearing in Eqs. (14) and (16) is determined from Darcy's law [52,64,63,66]: where k is the coefficient of permeability of the porous medium, g is the gravitational acceleration, and the density q fm is given by Eq. (7). For reasons of simplicity, the permeability in Eq. (19) is assumed to be isotropic. Note, however, that an extension towards an anisotropic permeability behaviour can be made in a straightforward fashion, and does not affect the basic features and structure of the model formulation. In accordance with the solid mechanics sign convention, the pore pressure p fm is assumed to take negative values, as a result of which the term in the right-hand side of Eq. (19) has a positive sign, i.e., the direction of the flux is in correspondence with the direction of an increasing (less negative) pore pressure. Since the porosity / fm of the granular medium increases during the erosion process, the coefficient of permeability k also increases, which is accounted for by means of a power law expression [54] k where k 0 is the initial permeability of the uneroded granular material, / fm;0 is the initial porosity, and the exponent m is a calibration factor. The value of m can be determined by matching Eq. (20) to the widely accepted Kozeny-Carman permeability relation [14,35,63] For example, for an initial porosity of / fm;0 ¼ 0:3, Eqs. (20) and (21) are in good correspondence over a large range of porosities 0 / fm 0:6 if the calibration parameter in Eq. (20) is set as m ¼ 4:6. This value is indeed selected for the numerical analyses presented in Sect. 5 and 6. The reason for using Eq. (20) instead of (21) in the present study is that Eq. (21) becomes singular when the soil material approaches a fully eroded state with / fm ! 1, which may induce numerical convergence problems.

Elasto-plastic constitutive model with the effect of erosion
The model for soil erosion defined in the previous section is now incorporated in a constitutive formulation for a cohesive-frictional granular material. The elasto-plastic deformation behaviour due to frictional sliding and granular compaction is modelled in accordance with the flow theory of plasticity, of which the basic concepts can be found in various text books, see e.g., [64]. After extending the elasto-plastic soil model with the effect of erosion, the isotropic hardening behaviour of the model is presented, followed by a discussion of the time integration scheme and the experimental calibration of the elastic and plastic material parameters.

Elasto-plasticity with erosion
Under the assumption of small strains, the total strain e of the granular skeleton may be additively decomposed into a reversible, elastic part e e and an irreversible, plastic part e p as During the erosion process the stiffness of the granular skeleton degrades, in accordance with the constitutive relation: with the erosion degradation parameter D er given by Eq. (18). The fourth-order isotropic elastic stiffness tensor D is determined by the Young's modulus E and Poisson's ratio m of the uneroded (D er ¼ 0) soil. Note from Eq. (18) that the uneroded soil is characterized by the initial particle volume fraction / s;0 . From the constitutive form given by Eq. (23), the degraded Young's modulus E er of the eroded material can be expressed as The stiffness reduction in accordance with Eq. (24) is analogous to a continuum damage mechanics formulation, in which the stiffness of a material with multiple small defects follows from the original stiffness of the undamaged material through a multiplication by ð1 À DÞ, with the damage parameter D reflecting the actual, total size of the defects [39].
In correspondence with the flow theory of plasticity, the plastic strain appearing in Eq. (23) follows from the time integration of the plastic strain rate as where _ j represents the magnitude of the plastic strain rate and m reflects the direction of plastic flow, Here, g is the plastic potential function, which will be defined further in this section. For the time integration of the effective stress, r 0 ¼ R t 0 _ r 0 ds, the plastic strain rate presented by Eq. (25) is substituted into the stress rate _ r 0 , which follows from Eq. (23) as The rate of degradation _ D er is hereby obtained from Eq. (18) as with _ / s determined by the mass-balance equation of the solid phase, Eq. (13) 1 . For the value b ¼ 1:0 selected in the present study, Eq. (28) simplifies to _ D er ¼ À _ / s =/ s;0 . Furthermore, the loading and unloading conditions of the elasto-plastic material with erosion are defined by the standard Kuhn-Tucker relations [36] where f represents the yield function. In accordance with the formulation proposed in [25], the yield function is composed of a frictional contour F fr that accounts for frictional sliding of particles, and a cap F c that characterizes volumetric compaction, see also Fig. 2: whereby the subindices ''fr'' and ''c'' of the frictional contour F fr and compression cap F c refer to ''friction'' and ''compaction'', respectively. Further, q 0 is the second deviatoric invariant of the effective stress, which is given by the usual definition in which s 0 is the deviatoric effective stress in the particle skeleton and p 0 is the hydrostatic invariant of the effective stress, where ''tr'' denotes the trace of the second-order (stress) tensor. In Eq. (30) the function C ¼ĈðhÞ, with h 2 ½0; p=3 representing the Lode angle, accounts for the variation in strength in the deviatoric plane in the transition from triaxial extension (h ¼ 0) to triaxial compression (h ¼ p=3). This function is given by [27] C ¼ whereby the triaxiality a defines the ratio between the strengths in triaxial extension and triaxial compression. The Lode angle h appearing in Eq. (33) satisfies the relation [64] cos ð3hÞ ¼ 27J 0 where the second deviatoric stress invariant q 0 is given by Eq. (31) and the third deviatoric stress invariant J 0 3 reads with ''det'' the determinant of the second-order (stress) tensor, and the deviatoric stress s 0 as in Eq. (31). As illustrated in Fig. 2a, the frictional contour F fr appearing in Eq. (30) is defined by a Drucker-Prager cone that has been extended with a smooth transition towards the tensile regime, in accordance with [23] F fr ¼ À 6 sinðuÞ 3 À sinðuÞ where u and c are the friction angle and cohesion of the soil, respectively, and c 1 and c 2 are calibration parameters. The current study focuses on erosion processes in a noncohesive sandy material, for which c ¼ 0. Nonetheless, for cohesive granular materials, the effect of erosion on the cohesive strength may be accounted for in a similar fashion as shown in Eq. (24) for the Young's modulus, i.e., c er ¼ ð1 À D er Þc (or an adapted form of this expression). The degraded cohesion c er then replaces the cohesion c in the current frictional failure model, see also [46,53]. Although, in principle, a similar approach can be followed to account for the effect of erosion degradation on the friction angle u, for simplicity it is assumed here that the friction angle remains constant and that the loss in material resistance during erosion is fully captured by the stiffness degradation, Eq. (23). This assumption is reasonable, since the loss in mechanical resistance resulting from relating erosion degradation to the material stiffness is comparable to relating it to the friction angle; in particular, at a given deformation the stress in the granular material in both approaches is forced to monotonically decrease towards zero when erosion develops towards completion, see also [53]. The cap model F c simulating volumetric compaction is illustrated in Fig. 2b, and has the form [25] Fig. 2 The yield function f of the elasto-plastic soil model is composed of two parts; a the contour F fr representing frictional sliding, Eq. (36), is combined with b the cap F c representing volumetric compaction, Eq. (37), to obtain c the overall yield function f given by Eq. (30), see also [25] The value of F c ranges between 1 (for p 0 ! p c ) and 0 (for p 0 ¼ X), see Fig. 2b, whereby at p 0 ! p c compaction does not play a role and the yield function, Eq. (30), becomes purely governed by the frictional surface, while at p 0 ¼ X the cap bounds the yield contour along the hydrostatic axis. A smooth connection between the cap and the frictional surface at the pressure value p 0 ¼ p c can be warranted by relating the frictional surface F fr;0 ¼F fr ðu ¼ u 0 ; p 0 ¼ p c Þ (i.e., the frictional surface evaluated at the initial friction angle u ¼ u 0 and the effective hydrostatic stress p 0 ¼ p c ) and p c to X as [25] where R is a scaling parameter. Note that the value of X develops because of the hardening behaviour of the cap, as described by the evolving compaction threshold p c . The specific evolution of p c is defined in Sect. 3.2 below. The flow rule presented by Eq. (26) requires the definition of a plastic potential g for determining the direction of plastic deformation. The plastic potential is taken as non-associative in order to adequately account for the plastic volumetric deformations (compaction and dilation) in the soil [64], whereby the specific form is assumed to be similar to that of the yield function, Eq. (30), but with the friction angle u in the frictional surface given by Eq. (36) replaced by the dilatancy angle w. Further, the Lode angle dependency C is ignored in the plastic potential by setting this contribution equal to unity, so that the plastic potential becomes characterized by a Drucker-Prager cone. This assumption is motivated from the experimental results reported in [33], which indicate that for frictional materials the direction of plastic deformation in the deviatoric plane corresponds reasonably well with a flow rule based on a Drucker-Prager cone, especially at moderate to low stress levels. Accordingly, the plastic potential has the form in which

Isotropic hardening
During frictional hardening, the isotropic evolution of the frictional surface F fr given by Eq. (36) is considered to be governed by the internal state variable j fr . Correspondingly, the friction angle u characterizing the frictional surface is assumed to depend on j fr by means of an exponentially-decaying law [56] u with u 0 the initial friction angle, u m the final, maximum friction angle and f fr a parameter defining the hardening rate. In a similar fashion, the dilatancy angle w appearing in the plastic potential g given by Eqs. (39) and (40) develops with j fr as [56] w where w 0 is the initial dilatancy angle and w m is the maximum dilatancy angle. In addition, the isotropic hardening response of the compaction cap F c defined by Eq. (37) is governed by the internal state variable j c . Correspondingly, the compaction threshold p c characterizing the compression cap evolves as where p c 0 is the initial value of p c and f c is the compaction rate. The effect of the internal state variables j fr and j c on the generated plastic deformation e p follows from decomposing the plastic flow rule, Eq. (25), into a frictional contribution and a compaction contribution where m fr and m c are the plastic flow directions related to frictional sliding and volumetric compaction, respectively. The frictional flow direction m fr can be obtained from Eqs. (26) and (39) through switching off the compaction contribution F c in Eq. (39) by equating it to unity, i.e., Similarly, the compaction flow direction m c is calculated from Eqs. (26) and (39) by switching off the frictional contributions in Eq. (39), which gives In order to relate _ j fr and _ j c to _ j, the plastic strain rate _ e p is decomposed into a deviatoric contribution _ c p and a volu- Combining Eq. (47) and the stress decomposition in Eq. (31) with Eqs. (25) and (26) allows to write _ c p and _ e p vol as In the same fashion, combining Eqs. (47) and (31) with Eqs. (44), (45) and (46) gives for _ c p and _ e p vol : Equating the deviatoric parts, Eqs. (48) 1 and (49) 1 , leads to Using this result in equating the volumetric parts, Eqs. (48) 2 and (49) 2 , gives with g, g fr and g c presented in Eqs. (39), (45) and (46), respectively. Eqs. (50) and (51) indicate how the rates of the internal state variables j fr and j c are related to the rate of j, as a result of which the evolution of the combined frictional-compaction yield function, Eq. (30), under isotropic hardening can be monitored via a single internal state variable j.

Time integration of the elasto-plastic erosion model
The time integration of the above elasto-plastic erosion model is performed using an incremental-iterative update procedure based on a fully implicit Backward Euler algorithm, with a consistent tangent operator formulated to construct the stiffness matrix at system level. In order to avoid the computation of complex, analytical derivatives, the consistent tangent operator is computed numerically using a perturbation method [49,59]. The numerical integration scheme has a similar structure as described in [25,56,59], i.e., the governing model equations are assembled in a vector of residuals, whereby the corresponding essential variables are solved iteratively for each time increment using a Newton-Raphson procedure. In this way, in an arbitrary material point of the simulated domain, the balance of linear momentum with b the body force per unit volume (e.g., as generated under dead weight loading), is satisfied under the appropriate boundary conditions.

Calibration of elasto-plastic material properties
The elastic and plastic material parameters of the soil model were calibrated using experimental data obtained from triaxial compression tests on a dense, well graded sand with some gravel, as reported in [58]. For this purpose, the uniform stress state generated in the triaxial compression test was mimicked by means of an FEM model of a single axisymmetric element, whereby the response of the elasto-plastic constitutive model was computed in accordance with the incremental-iterative update procedure described in Sect. 3.3. The material parameters of the elasto-plastic model were calibrated via a trial-and-error procedure, with the values listed in Table 1 leading to a good correspondence with the experimental data.
It is noted that in the triaxal compression test results presented in [58] the elastic stiffness characterizing the initial sample response is slightly dependent on the applied confining pressure; accordingly, the specific value of the Young's modulus presented in Table 1 on average matches the initial sample responses obtained in the range of confining pressures considered in [58]. Further, it can be concluded from Eqs. (41) and (42) that the calibrated values listed for the minimum and maximum dilatancy and friction angles result in w\u, which confirms that the model meets the second law of thermodynamics [64]. The maximum hydrostatic stress applied in the triaxial tests in [58] was p ¼ À157 kPa, at which the sample failed by frictional sliding without noticeable particle breakage. Hence, the values for the parameters p c 0 and R defining the compaction cap were determined based on additional experimental observations that particle crushing of noncohesive granular materials under one-dimensional compression is initiated at a hydrostatic pressure of about X ¼ À350 kPa [28,58]. In accordance with Eq. (38), at the onset of frictional sliding (whereby u ¼ u 0 ) this value can be matched when p c 0 ¼ À163 kPa and R ¼ 0:95, see Table 1. Note that the magnitude of the calibrated value of p c 0 ¼ À163 kPa is indeed larger than that of the maximum hydrostatic pressure of À157 kPa applied in the triaxial compression tests in [58]. Due to a lack of experimental data, the value of the compaction hardening rate was taken the same as that of the frictional hardening rate, f c ¼ f fr ¼ 250. Finally, the triaxiality a was determined from DEM simulation results of true triaxial tests performed on a cuboidal granular sample composed of equisized spherical particles, as presented in [57]. For a sample composed of freely rotating particles that have a local contact friction angle of 24 o , the maximum effective friction angle equals u m ¼ 22 o and the triaxiality was computed in [57] as a ¼ 0:82, while for a sample characterized by fully constrained particle rotation, the maximum effective friction angle is u m ¼ 50 o and the triaxiality was calculated as a ¼ 0:73. Using a linear interpolation, from the above values the triaxiality related to the present maximum friction angle of u m ¼ 40 o then becomes a ¼ 0:76, see Table 1.
With the above procedure, all elasto-plastic material properties are calibrated, by which the constitutive behaviour of the sand material is completely defined. For the practical case studies presented in Sect. 6, the influence of soil plasticity on the erosion profiles and surface settlements computed is identified by means of a comparison with numerical results obtained by using a simplified, elastic erosion model. The elastic erosion model is straightforwardly obtained from the above elasto-plastic formulation by setting the magnitude of the ''yield strength'' -i.e., the cohesion c in the frictional yield contour, see Eq. (36), and the pressure p c 0 required for initiating the compression cap, see Eq. (43), -artificially high.

Pore pressure field and staggered solution procedure
The coupled hydro-mechanical erosion problems presented in Sects. 5 and 6 are solved in an incremental-iterative fashion using a staggered numerical update scheme. In this section the pore pressure field is formulated from the storage equation derived in Sect. 2, and subsequently included in the staggered numerical update scheme. After that, the numerical implementation of the staggered update scheme is presented.

Pore pressure field
In order to derive the governing differential equation for the pore pressure field of the fluidic mixture, the storage equation given by Eq. (14) is used as a start: In correspondence with the solid mechanics sign convention adopted, see Sect. 2.1, in Eq. (53) a positive net outflow of the pore fluid ($ Á q fm [ 0) corresponds to a volume decrease ($ Á v\0), and a negative net outflow ($ Á q fm \0) corresponds to a volume increase ($ Á v [ 0). Using the strain decomposition in Eq. (22), the total volumetric strain rate in Eq. (53) may be decomposed into an elastic part and a plastic part as with the plastic volumetric strain rate _ e p vol given by Eq. (48) 2 . The elastic volumetric strain rate _ e e vol can be computed by first specifying the general constitutive expression, Eq. (23), in terms of the hydrostatic effective stress p 0 , Eq. (32), as where B is the elastic bulk modulus. Eq. (55) can be expressed in rate form as Further, from Terzaghi's stress decomposition, Eq. (17), the hydrostatic effective stress rate _ p 0 appearing in the lefthand side of Eq. (56) reads with the time rate change of the hydrostatic total stress as From Eqs. (56) and (57), the elastic volumetric strain rate is obtained as Combining Eqs. (54) and (59), and inserting the result, together with Darcy's law, Eq. (19), into Eq. (53), the nonlinear partial differential equation describing the pore pressure field in a poro-elasto-plastic medium with erosion becomes For a poro-elastic medium with erosion, the plastic deformation term in Eq. (60) vanishes ( _ e p vol ¼ 0), leading to For a poro-elastic medium without erosion with k 0 the initial permeability coefficient that appears in Eq. (20). When the transport of the fluidic mixture in the pores occurs under a constant total stress, as prescribed through fixed traction boundary conditions, the time rate of change of the total stress is zero, _ r ¼ 0, by which _ p ¼ 0, so that the right-hand side of Eq. (62) vanishes. Conversely, for poro-elastic systems characterized by displacement boundary conditions in one or more directions, the hydrostatic stress rate typically does not vanish, _ p 6 ¼ 0, which thus is accounted for via the right-hand side of Eq. (62). Note further that Eq. (62) describes the consolidation process of a poro-elastic medium, during which the two-phase material gradually changes volume in response to a change in pressure.
Consider now the specific case of one-dimensional consolidation of a poro-elastic medium under a constant total stress r in the (axial) x-direction. Since the one-dimensional medium does not deform in the two lateral directions, the hydrostatic effective stress p 0 ¼ p À p fm can be expressed in terms of the effective stress r 0 in axial direction as p 0 ¼ r 0 ð1 þ mÞ=ð3 À 3mÞ, where r 0 ¼ r À p fm . With these expressions, the hydrostatic stress rate is obtained from Eq. (57) as _ p ¼ _ p fm ð2 À 4mÞ=ð3 À 3mÞ. Inserting this result into the right-hand side of Eq. (62) turns this equation into Terzaghi's one-dimensional consolidation equation [52,62,66]: with c v representing the consolidation coefficient and m v the one-dimensional compressibility coefficient. In order to identify some important features and characteristic time scales of the soil consolidation and erosion processes, in Sect. 5 an analysis of a one-dimensional poro-elastic benchmark problem is performed in two separate steps. In the first step, a poro-elastic sample is subjected to a specific change in pore pressure and is allowed to consolidate without erosion. The FEM solution for this problem will be validated to the corresponding analytical solution of Terzaghi's one-dimensional consolidation equation, Eq. (63). In the second step, the one-dimensional flow problem is extended by including the effect of erosion. Accordingly, the pore pressure field is described by Eq. (61), with the permeability coefficient k in Eq. (61) being dependent on the actual porosity, see Eq. (20). For reasons of simplicity, the influence of material plasticity is left out of consideration in the analyses by adopting an artificially high yield strength in the constitutive soil model; however, as already mentioned, plasticity effects will be accounted for in the practical case studies analysed in Sect. 6, which relate to a concrete sewer pipe embedded in a sandy soil structure.

Staggered solution procedure
The coupled hydro-mechanical analyses were performed with the commercial finite element program ABAQUS Standard. 1 The problem is fully described by which makes the spatial discretization of the solid volume fraction / s unnecessary, see also [38]. The constitutive elasto-plastic erosion model presented in Sect. 3 has been implemented in ABAQUS as a user-supplied subroutine, i.e., a ''UMAT'', in accordance with the time integration scheme outlined in Sect. 3.3. Since Eq. (60) is analogous to the partial differential equation for the process of non-linear heat conduction, the thermalmechanical module in ABAQUS was used for carrying out the coupled hydro-mechanical analyses. By relating the material parameters of the thermo-mechanical model and hydro-mechanical model as derived in Appendix A, the temperature and thermal flux defining the thermo-mechanical process can be respectively interpreted as the pore fluid pressure and the pore fluid flux characterizing the hydro-mechanical process. The coupled hydro-mechanical response was solved by using the staggered incrementaliterative update scheme summarized in Table 2.
The hydrological and mechanical field variables were solved for each time increment in a sequential fashion, with the couplings between the individual fields accounted for through a temporal extrapolation. The time increment was chosen relatively small, such that the error introduced by the time discretisation had a negligible influence on the numerical result. Each time increment started with the computation of the hydrological field variables p fm and q fm from the mass balance equation for the fluidic mixture phase, Eq.   the volumetric plastic strain rate _ e p vol was obtained from Eq. (48) 2 , the degradation parameter D er and its rate _ D er were calculated from Eqs. (18) and (28), respectively, the volumetric elastic strain e e vol was determined via the inverse of Eq. (55), and the rate of eroded mass term _ m er was computed from Eq. (16). The pore pressure values and flux values were subsequently transferred to the mechanical analysis, after which the balance of linear momentum, Eq. (52), was solved together with the mass balance equation for the solid phase, Eq. (13) 1 , using the constitutive soil model and the iterative solution procedure as described in Sect. 3.3. For this purpose, the effective stress r 0 employed in the constitutive soil model was calculated by substituting the updated pore pressure field p fm in the stress decomposition, Eq. (17). From the converged solution, the displacement field u and the erosion and plasticity internal state variables / s and j were updated. The updated solid volume fraction / s was used to update the porosity / fm via Eq. (3), which was subsequently transferred to the hydrological analysis in order to update the permeability k through Eq. (20), as required for solving the pore pressure field p fm in the next increment via Eq. (60), and the flux field q fm via Eq. (19). Similarly, the updated values of the time rate of change of the total hydrostatic stress _ p (in accordance with Eq. (58)), the volumetric plastic strain rate _ e p vol (in accordance with Eq. (48) 2 ), the degradation parameter D er and its rate _ D er (in accordance with Eqs. (18) and (28), respectively), the volumetric elastic strain e e vol (in accordance with the inverse of Eq. (55)), and the rate of eroded mass _ m er (in accordance with Eq. (16)) were transferred to the hydrological analysis in order to solve the mass-balance equation, Eq. (60), in the next time increment. Subsequently, the above procedure was repeated for the next time increment.

One-dimensional benchmark problem
In order to demonstrate the basic features of the erosion model and validate the numerical implementation of the coupled hydro-mechanical formulation, in this section a one-dimensional poro-elastic benchmark problem is analysed in two steps. Firstly, a poro-elastic soil sample is subjected to an instantaneous change in pore pressure, whereby only the consolidation behaviour of the sample is studied through suppressing the effect of erosion. The FEM results are validated via a comparison with the analytical solution of the problem. Secondly, the poro-elastic flow problem is extended with the effect of erosion, which illustrates the specific influence on the degradation and deformation behaviour of the sample, and indicates the relation between the characteristic time scales for the erosion and consolidation processes. The model parameters will be presented first, followed by a discussion of the computational results of the benchmark problem.

Model parameters
The geometry and boundary conditions of the one-dimensional benchmark problem of a poro-elastic sample with length l ¼ 1 m are illustrated in Fig. 3. In the 2D FEM representation, the height of the sample is taken as h ¼ 0:01 m, whereby the one-dimensional character of the hydro-mechanical problem is mimicked by modelling the out-of-plane direction as plane strain, and prescribing the displacement and flux of the fluidic mixture in the normal direction of the horizontal top and bottom boundaries to be zero, u Á n ¼ 0 and q fm Á n ¼ 0, with n the unit outward normal vector at a boundary. Under these lateral boundary conditions, in the FEM simulations the lateral total stress does not remain constant during the simulation, as a result of which the rate of the total hydrostatic pressure is unequal to zero, _ p 6 ¼ 0. Accordingly, the FEM simulation needs to be performed by considering the differential equation, Eq. (62), for the poro-elastic problem without erosion, and by considering Eq. (61) for the poro-elastic problem with erosion.
The initial state of equilibrium of the poro-elastic sample (i.e., the equilibrium state at time t ¼ 0) is characterized by a uniform pore pressure of p fm;0 ¼ À20 kPa, which is imposed as an initial pore pressure applied to the sample. The initial axial total stress in the sample is also uniform, in Fig. 3 Two-dimensional, plane-strain finite element model of the onedimensional benchmark problem, with displacement constraints at the top, bottom and right boundaries of the domain, an initial traction boundary condition r 0 ¼ À100 kPa at the left boundary, and a zero outward flux condition at the top and bottom boundaries; at t [ 0 a horizontal pore water flow is initiated by abruptly increasing the magnitude of the pore pressure at the left boundary to p l fm ¼ À50 kPa while keeping the pore pressure p r fm at the right boundary equal to the initial pore pressure in the domain, p fm;0 ¼ À20 kPa accordance with the traction boundary condition, r 0 ¼ À100 kPa, applied at the left sample side, and a zero axial displacement, u r ¼ 0, at the right sample side. At t [ 0 the pore pressure at the left side of the sample is increased in magnitude towards a value of p l fm ¼ À50 kPa, while at the right side of the sample the initial value of the pore pressure is maintained, p r fm ¼ À20 kPa. Here, the superindices l and r respectively designate ''left'' and ''right''. This abrupt change in pore pressure is supposed to qualitatively mimic the effect of a local increase in groundwater level, e.g., as a result of heavy rainfall or a flood. Due to the local increase in pore pressure, an axial flux q fm is initiated at the left sample boundary, as a result of which the sample will expand in volume over time.
In the FEM model, the mechanical and hydrological responses are computed by means of coupled elements that both have temperature -which may be interpreted as the pore pressure, see Appendix A -and displacement degrees of freedom. The mesh is generated by 255 plane-strain 4node iso-parametric elements, equipped with a 2 Â 2 Gauss quadrature. The mesh is refined at the left side of the domain in order to accurately describe the relatively large spatial gradient in pore pressure that will occur under the abrupt pore pressure change imposed by p l fm , see Fig. 3b. The material parameters related to pore water flow and erosion are listed in Table 3, and are assumed to be representative of a dense, well graded sand with some gravel, i.e., a similar granular material as considered for the calibration of the elasto-plastic material properties in Section 3.4. Hence, the elastic parameters E and m of the (uneroded) sand material are taken from Table 1. Further, the values selected for the initial permeability k 0 and the initial volume fraction of solid particles / s;0 are adopted from [66]. The value of the parameter m is based on a calibration of the Kozeny-Carman relation, Eq. (21), as explained in Sect. 2.1. Due to the absence of experimental data, the parameters of the erosion model are estimated from engineering judgement. For the present benchmark problem it is assumed that erosion takes place in accordance with a suffusion process. Accordingly, only finer soil particles are eroded, whereby the particle contact force fabric, and thus the bearing strength of the sample, to some extent is maintained, in accordance with the selection of a final, minimum value of the solid volume fraction of / min s ¼ 0:35. The value / min s ¼ 0:05 that is representative of a suffosion type of erosion process, which is also listed in table 3, will be used in a practical case study presented in Sect. 6. Together with an initial solid volume fraction of / s;0 ¼ 0:7, the value / min s ¼ 0:35 leads to a maximum degradation parameter of D max er ¼ ð0:70 À 0:35Þ=0:70 ¼ 0:5, see Eq. (18). Hence, the final elastic stiffness E er of the degraded soil material is ðE er =EÞ Â 100% ¼ 50% of the initial elastic stiffness E, see Eq. (24). The maximum permeability reached at D max er ¼ 0:5 equals k max ¼ 3:5 Â 10 À3 m/s, see Eq. (20), which is indeed representative of a relatively loose sand [66]. For long-term sand erosion processes -in the order of days or (much) longer -the concentration of fluidized solid particles in the fluidic mixture typically is very low, c fs ( 1; accordingly, this value has been selected as c fs ¼ 0:01, which, with Eq. (7), results in a fluidic mixture density of q fm ¼ 1014 kg/ m 3 . Hence, the fluidic mixture density q fm is only 1:6% larger than the pore water density q f ¼ 997 kg/m 3 , so that the results computed for the pore pressure p fm and flux q fm may be essentially interpreted as ''pore water pressure'' and ''pore water flux'' fields, respectively. Finally, the gravitational acceleration equals g ¼ 10 m/s 2 .

One-dimensional consolidation without erosion
The FEM results calculated for one-dimensional consolidation without erosion are compared with the corresponding analytical solution. The analytical and numerical results are evaluated and compared by considering the time evolution of 4 field parameters, namely the pore pressure, p fm ¼p fm ðx; tÞ, the axial flux, q fm ¼q fm ðx; tÞ, the axial displacement, u ¼ûðx; tÞ, and the axial effective stress, r 0 ¼r 0 ðx; tÞ. The analytical solution for the pore pressure field p fm ¼p fm ðx; tÞ is obtained by solving Eq. (63) under the specific initial pressure field p fm;0 and pressure boundary conditions p l fm and p r fm applied, see also Fig. 3. The closed-form expression for p fm ¼p fm ðx; tÞ has been derived following a standard solution procedure for non-homogeneous, linear partial differential equations [52]. Substituting this solution into Eq. (19) results in the axial flux field q fm ¼q fm ðx; tÞ. In accordance with Terzaghi's stress decomposition, Eq. (17), the axial effective stress is subsequently calculated as r 0 ¼r 0 ðx; tÞ ¼ r 0 Àp fm ðx; tÞ. Finally, the axial displacement field is obtained from the spatial integration of the mechanical constitutive relation, u ¼ûðx; tÞ ¼ m v Rr 0 ðx; tÞ À r 0 0 dx, and is thus evaluated by taking the sample with the initial stress, r 0 0 ¼ r 0 À p fm;0 , as the zero reference state. The integration constant following from the integration procedure is set by the displacement boundary condition u r applied at the right sample side. Accordingly, the analytical solution of the pore pressure field p fm ¼p fm ðx; tÞ is obtained as and the axial flux q fm ¼q fm ðx; tÞ has the form In addition, the axial effective stress field r 0 ¼r 0 ðx; tÞ follows as and the axial displacement field u ¼ûðx; tÞ reads with C 1 ¼Ĉ 1 ðtÞ given by In the above expressions, the consolidation coefficient c v and the compressibility coefficient m v are given by Eq. (63). The total number of terms used for constructing the summation series in Eqs. (64) to (68) is selected as n ¼ 20, which turned out to be sufficient for obtaining a converged analytical result. Figures 4a,b,c, and d illustrate the analytical and numerical solutions for, respectively, p fm , q fm , r 0 and u, as considered across the sample length l ¼ 1:0 m at various time instants. It can be observed that the analytical and numerical solutions are in perfect agreement. As a result of the applied boundary conditions, the pore pressure p fm in the sample gradually evolves from the uniform reference value p fm ¼ p fm;0 ¼ À20 kPa at t ¼ 0 to a linear, steadystate profile at t ¼ 2:0 s, see Fig. 4a. The steady-state pore pressure field p fm;ss ¼p fm;ss ðxÞ is described by the first two terms in the right-hand side of Eq. (64), i.e., and is reached relatively fast, namely at t ¼ 2:0 s, which is due to the high permeability k 0 ¼ 10 À4 m/s of the well graded sand material. The profile of the flux depicted in Fig. 4b develops non-uniformly across the sample, with the value initially raising from zero at the right sample boundary towards a maximum at the left sample boundary (which is where the pore pressure is maximal). Hence, the fluidic mixture flows inwardly from the left sample boundary, whereby the local peak value of the flux gradually decreases with time, eventually leading to a uniform flux profile at steady state. The expression for the flux q fm;ss reached at steady state, t ¼ 2:0 s, follows from the first term in Eq. (65) as which, as can be confirmed from Fig. 4b, corresponds to a value q fm;ss ¼ 3:0 Â 10 À4 m/s. Figure 4c shows that the axial effective stress, which is initially uniform in accordance with the initial condition, r 0 0 ¼ r 0 À p fm;0 ¼ À100 þ 20 ¼ À80 kPa, grows non-uniformly as a result of the different left and right boundary values, i.e., r 0 l ¼ r 0 À p l fm ¼ À100 þ 50 ¼ À50 kPa and r 0 r ¼ r 0 À p r fm ¼ À100 þ 20 ¼ À80 kPa. The spatial variation of the axial effective stress asymptotes with time towards a linear stress profile at steady state, r 0 ss ¼r 0 ss ðxÞ, which is given by the first three terms in the right-hand side of Eq. (66): Clearly, from the elastic constitutive relation the linear steady-state stress profile translates into a linear steadystate strain profile, and thus into the quadratic steady-state displacement profile depicted in Fig. 4d for t ¼ 2:0 s. The expression for the steady-state displacement field u ss 1 u ss ðxÞ is obtained from Eqs. (67) and (68) as Using the specific material parameters and boundary conditions, with Eq. (72) the steady-state displacement at the left sample boundary is calculated as u ss ð0Þ ¼ À0:93 mm. Note that this negative amplitude can be confirmed from the response at t ¼ 2:0 s as depicted in Fig. 4d. The time development of the displacement pattern from a uniform zero value at t ¼ 0 to a quadratic steady-state displacement profile with a negative amplitude at t ¼ 2:0 s illustrates that the sample swells under the instantaneous increase in pore pressure applied at its left boundary. Fig. 3 for the case of an elastic sandy soil that consolidates without erosion; for a selection of time instants (measured in seconds), the figure shows the spatial variation of a the pore pressure p fm , b the axial pore water flux q fm , c the axial effective stress r' and d the axial displacement u

One-dimensional consolidation with erosion
The effect of erosion in the one-dimensional flow problem sketched in Fig. 3 is first analysed by considering the time evolutions of the rate of the eroded mass _ m er (scaled by q s ) and the degradation parameter D er , see Figs. 5a and b.
It can be observed that erosion becomes noticeable after the consolidation process has reached a steady state at t ¼ 2:0 s. Since the rate of eroded mass _ m er is proportional to the flux q fm , see Eq. (16), and the flux at steady-state consolidation is more or less uniform across the sample, see Fig. 6b (hereby disregarding the minor spatial gradient that quickly vanishes for t [ 2:0 s), the erosion process develops spatially uniformly in time. Figure 5b shows that the degradation parameter D er evolves relatively slowly, which indicates that the characteristic time scale of the   Fig. 3 for the case of an elastic sandy soil that consolidates with erosion; for a selection of time instants (measured in seconds), the figure shows the spatial variation of a the pore pressure p fm , b the axial pore water flux q fm , c the axial effective stress r 0 and d the axial displacement u erosion process is much larger than that of the consolidation process. Note further that at t ¼ 5 Â 10 5 s (= 5.8 days) the suffusion type of erosion has become maximal, in correspondence with a maximum degradation parameter of D max er ¼ 0:5. Hence, the erosion process has finished, whereby the rate of eroded mass _ m er in Fig. 5a has uniformly dropped to zero, in line with the second expression in Eq. (16). Figures 6a,b,c, and d respectively show the pore pressure p fm , the axial flux q fm , the axial effective stress r 0 and the axial displacement u during the erosion process.
As the consolidation process is in a steady state, the linear pore pressure profile p fm given by Eq. (69) does not change during erosion, see Fig. 6a. Conversely, Fig. 6b illustrates that the flux q fm becomes larger with increasing erosion, which obviously is due to the increasing permeability of the sand material, see Eqs. (19) and (20). The maximum flux reached at the end of the erosion process, t ¼ 5 Â 10 5 s, can be calculated from Eqs. (19) and (20) as q max fm ¼ k max ðp r fm À p l fm Þ=ðq fm glÞ ¼ 1:04 Â 10 À2 m/s, which is in agreement with the value depicted in Fig. 6b. Further, in accordance with Eq. (16), the increase of the flux q fm generates an increase of the rate of eroded mass _ m er , as depicted in Fig. 5a. The steady-state effective stress profile r 0 shown in Fig. 6c obviously remains constant during the erosion process, in accordance with Terzaghi's stress decomposition, Eq. (17). The displacement profile shown in Fig. 6d monotonically grows during the erosion process. The final profile at the end of the erosion process, t ¼ 5 Â 10 5 s, is characterized by a spatially uniform degradation parameter of D er ¼ 0:5 that effectively reduces the sample stiffness to E er ¼ 0:5E, see Eq. (24). Hence, the analytical form of the displacement profile at t ¼ 5 Â 10 5 s follows from multiplying the steady-state displacement profile at full consolidation, Eq. (72), by a factor of E=E er ¼ 2, which leads to a maximum displacement value at the left sample boundary of 2u ss ð0Þ ¼ À1:86 mm, see Fig. 6d. In other words, the erosion process causes that the swelling of the consolidated sample increases by a factor of two.
6 Practical case studies on soil erosion near a sewer system In this section two practical case studies will be considered, which relate to a sewer system embedded in a sandy soil structure. The first case study simulates the process of soil piping caused by suffusion, as occurring near a sewer system subjected to natural groundwater flow. The second case study models void formation caused by suffosion, as generated under a strong groundwater flow towards a defect sewer pipe, i.e., a sewer pipe with a gap created by a sudden failure of a pipe connection. The model parameters will be presented first, followed by a discussion of the computational results of the two case studies.

Model parameters
The specific geometry considered in the hydro-mechanical analyses is depicted in Fig. 7.The soil structure consists of a top layer of dry sand with a thickness of 2.27 m, which is supported by a fully saturated sand layer with a thickness of 1.50 m, followed by a relatively thick, fully saturated clay layer. The location of the centre of the sewer pipe corresponds with the ground water table at 2.27 m depth. The inner diameter and wall thickness of the round sewer pipe are 400 mm and 65 mm, respectively. The sewer system is covered by 2 meters of dry sand, which is representative of the conditions in the Netherlands [9]. The specific part of the geometry simulated with FEM is indicated in Fig. 7 by the light shaded rectangular section of 6 Â 3 m 2 . The mechanical and hydrological initial and boundary conditions used in the FEM simulation are  Fig. 8a for the first case study, and in Fig. 8b for the second case study. The finite element discretization of the simulated domain is shown in Fig. 8c. The ground water table initially equals the phreatic surface at which the pore pressure is zero, p fm ¼ 0. The hydrological boundary condition imposed at the ground water table corresponds to a zero normal flux q fm Á n ¼ 0, so that the ground water does not migrate into the dry sand layer above. With this boundary condition, the phreatic surface during the simulation is allowed to move in the downward direction, thereby creating a capillary area in between the ground water table and the phreatic surface in which matric suction takes place, i.e., the pore water pressure becomes positive (in accordance with the solid mechanics sign convention adopted in this study). This situation, which has some relevance in the second practical case study, has been checked a priori, showing that the capillary zone characterized by matric suction typically remained relatively small, i.e., the distance between the ground water table and the phreatic surface always appeared to be less than one half of the pipe diameter, whereby the maximum value of matric suction was ?2.4 kPa, which is a realistic value for a sandy soil [26,72].
The clay layer located below the sand layer has a coefficient of permeability that is typically two to four orders of magnitude lower than that of the sand layer [66], so that it may be assumed as impermeable. Accordingly, at the bottom of the saturated sand layer a zero normal flux boundary condition is adopted, q fm Á n ¼ 0. In a similar fashion, the outer circumference of the concrete sewer pipe is modelled as impermeable.
In the first case study, at t [ 0 the horizontal groundwater flux normal to the left domain boundary is set equal to the erosion threshold value, q fm Á n ¼ q fm;cr ¼ 10 À6 , see also Table 3. At the right domain boundary the dynamic pore pressure is set to zero, so that the total pore pressure is Fig. 8 FEM models of the two practical case studies, with a the mechanical/hydrological initial and boundary conditions of the first practical case study representative of a sewer system subjected to natural ground water flow, b the mechanical/hydrological initial and boundary conditions of the second practical case study representative of a strong ground water flow near a defect sewer pipe, and c the finite element discretization used in the two case studies; the vertical traction r 0 along the upper domain boundary reflects the dead weight of the 0.77 m of dry sand material located above the computational domain, see Fig. 7 maintained at the geostatic pore pressure, p fm ¼ p fm;0 , which increases linearly with depth. From Eqs. (15) and (19), it may be concluded that these boundary conditions relate to steady-state values for the groundwater velocity and the pressure gradient that are realistic for typical in-situ conditions [20,45,55]. Under the application of the above boundary flux a natural groundwater flow is initiated from the left to the right domain boundary, whereby at locations near the circular sewer pipe the flux is expected to increase, thereby exceeding the erosion threshold value, kq fm k [ q fm;cr .
In the second case study, at t [ 0 the sewer pipe experiences at its bottom a defect with a width of 20 mm, i.e., a gap created by a sudden failure of a pipe connection. Across the width of the gap, the pore pressure is set equal to the inner pressure in the pipe, which agrees with the atmospheric pressure, p fm ¼ 0 kPa. The gap width corresponds to the maximum grain size present in the grain-size distribution of the sandy material (with some gravel) tested in [58], which, in accordance with the parameter values in Table 1, characterizes the dry and saturated sandy layers of the soil configuration depicted in Fig. 7. Hence, all particles of the sand material in principle can flow through the gap opening. At the left and right domain boundaries the pore pressure is maintained at the initial geostatic pore pressure, p fm ¼ p fm;0 . Accordingly, a relatively strong groundwater flow is initiated from the soil structure into the defect sewer pipe, with the chance of locally washing away soil material and creating a suffosion erosion void.
In both case studies, the initial stresses rðt ¼ 0Þ ¼ r 0 in the simulated domain are generated by the dead weight loading of the soil layers and the sewer system. Accordingly, in each material point of the modelled soil structure the geostatic effective initial stress r 0 0 is computed from the total initial stress and the geostatic pore pressure in accordance with Eq. (17), r 0 0 ¼ r 0 À p fm;0 1. The upper boundary of the computational domain is subjected to a vertical normal stress r 0 that reflects the dead weight of the 0.77 m of dry sand above. At the lower boundary of the domain, which represents the bottom of the saturated sand layer, the vertical displacement is prescribed as zero, and at the left and right domain boundaries the horizontal displacement is set to zero. With these boundary conditions, the initial stresses r 0 in the domain are calculated by applying the dead weight loading incrementally in a sequence of small equilibrium steps. Table 4 summarizes the initial dry and wet densities of the sand. For the dry sand material, the initial densitỹ q dry;0 ¼ 1855 kg/m 3 was computed in accordance with Eq. (5), by multiplying the density of quartz, q s ¼ 2650 kg/ m 3 , with the initial volume fraction of particles, / s;0 ¼ 0:7 listed in Table 3. For the wet, saturated sand this density was increased with the partial density of the available pore water, ð1 À / s;0 Þq fm ¼ 299 kg/m 3 , leading toq wet;0 ¼ 2154 kg/m 3 . The material properties of the concrete sewer pipe were obtained from Eurocode EN1992-1-1, and correspond to a concrete strength class C50/60. The constitutive behaviour of the concrete pipe is simulated by means of a basic linear elastic model, which is considered to adequately reflect the mechanical interaction of the pipe with the surrounding soil. It is hereby mentioned that a detailed study of the stress redistributions in the pipe under the development of the soil erosion falls beyond the scope of the present work; more details on this aspect can be found in [51].
The elasto-plastic properties of the sand material are listed in Table 1, and the material parameters related to pore water flow and erosion in the sand are summarized in Table 3. Similar to the benchmark problem studied in Sect. 5, in the first case study the minimum volume fraction of sand particles is taken / min s ¼ 0:35. This value corresponds to a maximum degradation parameter of D max er ¼ ð/ s;0 À / min s Þ=/ s;0 ¼ 0:35=0:7 ¼ 0:50, and is considered to be representative of a suffusion type of erosion. Accordingly, some particle contact force fabric is maintained, such that the final, minimum elastic stiffness within the eroded soil pipe geometry is half the initial elastic stiffness of the sand, E er ¼ 0:5E, see Eq. (24). Conversely, in the second case study the minimum volume fraction of solid particles is set considerably smaller, in correspondence with the value listed in Table 3 for the case of suffosion, i.e., / min s ¼ 0:05, whereby the maximum degradation parameter corresponds to D max er ¼ 0:93. Hence, at this stage the soil structure supporting the defect sewer pipe will be almost completely washed away by the strong ground water flow, thereby creating an erosion void for which the minimum stiffness is close to zero, i.e., E er ¼ 0:07E. Setting the value of / min s for suffosion erosion slightly larger than zero thus prevents the minimum material stiffness from becoming zero, as a result of which the structural stiffness matrix of the FEM model remains well-conditioned during the numerical solution procedure, and convergence of the solution can be warranted.
In each loading step of the simulations, the wet density of the sand is updated in accordance with Eqs. (3) to (5) using the actual particle volume fraction / s of the sand. Hence, when the volume fraction of solid particles / s has reached its minimum value of / min s ¼ 0:35 and / min s ¼ 0:05 in, respectively, the first and second case study, the wet density becomes minimal, with the specific value for the ''degraded sand material'' computed from Eqs. (3) to (5) as q wet ¼ 1576 kg/m 3 andq wet ¼ 1080 kg/m 3 , respectively.
The computational domain is modelled as two-dimensional, assuming a plane-strain condition in the out-ofplane direction. The soil structure with the elasto-plastic sand behaviour is discretized with 47046 plane-strain 6-node iso-parametric coupled temperature-displacement elements, which are equipped with a 3-point Gauss quadrature. The sewer pipe is simulated as elastic, and is discretized by 2089 plane-strain 6-node iso-parametric elements with a 3-point Gauss quadrature. As indicated in Fig. 8c, the mesh of the soil structure is refined towards the (coherent) interface with the sewer pipe in order to accurately simulate the local erosion profile. A preliminary mesh refinement study has indicated that the present discretization is sufficiently fine for obtaining converged numerical results for the two case studies. 6.2 Soil piping due to natural groundwater flow near a sewer system The FEM simulation starts with the application of the initial loading generated by the dead weight of the soil structure and the concrete pipe. The vertical settlement profile resulting from the dead weight loading is depicted in Fig. 9a, and the corresponding deviatoric plastic deformations are illustrated in Fig. 9b. The plastic deformations are reflected in the contour plot by the deviatoric invariant c p ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi  Fig. 9a that the initial settlement profile is fairly uniform along the horizontal direction, with the settlements above the sewer pipe being slightly smaller than at the left and right domain boundaries. Obviously, this difference may be ascribed to the relatively low effective density of the hollow sewer pipe. As indicated in Fig. 9b by the grey-hatched area, the deviatoric plastic deformation is maximal at the top-left and top-right of the pipe circumference, reaching a value of c p ¼ 0:45%. It can be further observed that the whole domain undergoes plastic deformations after the application of the dead weight loading. The initial coefficient of lateral earth pressure for the non-cohesive sand can be computed from the ratio between the horizontal and vertical normal stresses, K ¼ r xx =r zz , and varies between between K ¼ 0:25 and K ¼ 0:36 going from the bottom to the top along the left (or right) boundary of the computational domain. In accordance with the range defined by the initial and maximum friction angles of the sand material as listed in Table 1, i.e., 29 o u 40 o , this means that the sand along the left (or right) domain boundary experiences a horizontal stress state that approaches the condition of active lateral earth pressure, as characterized by K ¼ ð1 À sin uÞ=ð1 þ sin uÞ [66]. In contrast, right next to the sewer pipe the vertical stress generated by the dead weight loading is considerably smaller than the horizontal stress, in correspondence with K ¼ 3:29, which is caused by a local arching of the sand material located along the upper half of the sewer pipe.
The application of a groundwater flux of q fm Á n ¼ 10 À6 m/s at t [ 0 at the left boundary of the computational domain initiates a groundwater flow from the left to the right domain boundary. Due to the relatively large permeability of the sand material, at sufficient distance below the sewer pipe the groundwater flow within a few minutes of time reaches an almost uniform steady-state profile in the horizontal direction, whereby the flux value is equal, or close to, the applied boundary flux kq fm k ¼ 10 À6 m/s, and the spatial pattern of kq fm k has become virtually symmetric with respect to the vertical center line of the FEM model. Note that a fast convergence towards a steady-state flow pattern was also observed for the one-dimensional benchmark problem analysed in Sect. 5.
The time evolution of the erosion profile developing under the groundwater flow is sketched in Fig. 10, by plotting the degradation parameter D er at three different time instants, namely 7 years (Fig. 10a), 13 years (Fig. 10b) and 20 years (Fig. 10c). It can be seen that the erosion   Fig. 13 Time evolution of the vertical displacement w (in mm) of the ground surface -measured above the centre of the sewer pipe -under natural groundwater flow (first case study) for i) an elasto-plastic sand material and ii) an elastic sand material profile starts to grow at the bottom of the sewer pipe, and subsequently spreads in more or less horizontal direction towards the left and right domain boundaries, thereby creating a soil piping profile. Note that the green colour in the contour plot designates the area in which the erosion degradation parameter has reached its maximum value, D er ¼ D max er ¼ 0:5, and that the characteristic time scale associated to the development of the erosion profile is much larger than that of the process of ground water flow, i.e., in the order of tens of years.
The groundwater flux profiles after 7, 13 and 20 years are depicted in Figs. 11a, b and c, respectively, using the Euclidian norm of the flux, kq fm k, as the contour plot variable. It can be confirmed that the area in which the groundwater flux kq fm k is larger than the threshold value of q fm;cr ¼ 10 À6 m/s indeed corresponds to the geometry of the erosion profile shown in Fig. 10. Further, the groundwater flux appears to be maximal directly below the pipe, whereby the flux value increases when erosion develops, as caused by the increase in permeability, see Eq. (20). Figure 12 depicts the contour plot with the deviatoric plastic deformations, as evaluated after the application of the dead weight loading and an erosion process of 20 years. For clarity, in the contour plot the same scale division is used as for the deviatoric plastic deformations generated after the application of only the dead weight loading, see Fig. 9b. In comparison with the deviatoric plastic deformations generated after the dead weight loading, left and right below the pipe the erosion process has led to a maximal increase in plastic strain of about a factor of two, in correspondence with a value of c p ¼ 0:74%. The evolution of the surface displacement, measured above the centre of the sewer pipe, is depicted in Fig. 13. In order to clearly identify the displacement contribution caused by plastic deformations, the surface response is compared to that from a simulation in which erosion takes place on a fully elastic soil material. As explained in Sect. 3.4, the elastic soil model is obtained from the elastoplastic formulation by setting the yield strength artificially high. The erosion profile following from the simulation with an elastic sand material is similar to that with an elasto-plastic sand material, see Fig. 10, from which it may be concluded that the contribution of the plastic volume change _ e p vol to the mass balance equation of the fluidic mixture, Eq. (60), is minor. Figure 13 illustrates that for an elasto-plastic sand material the surface displacement monotonically increases with the development of erosion, whereby the value reached after 20 years equals 0.60 mm. The response for an elastic soil material shows a similar trend, whereby the deformation reached after 20 years is 25% less, i.e., 0.45 mm. Since the surface displacements remain rather small, the practical suitability of using these as a monitoring parameter for the detection of a suffusion type of erosion under natural groundwater flow appears to be limited. This conclusion is in correspondence with other scientific studies, which report that soil pipes are only observable at the soil surface when a pipe roof collapses; accordingly, they are considered as ''apparently inactive'' over a long period of time, until clear surface evidence appears [2, 3, 65].

Void formation due to strong groundwater flow near a defect sewer system
Similar as in the first case study, the FEM analysis of soil erosion near the defect sewer pipe starts with the application of the dead weight loading of the soil structure and the sewer pipe. The atmospheric pressure p fm ¼ 0 kPa applied at t [ 0 as a boundary condition across the gap width at the bottom of the sewer pipe (see Fig. 8b) initiates a relatively strong groundwater flow into the sewer pipe. The erosion profile caused by this flow profile is depicted in Fig. 14 at three different time instants, namely 47 days (Fig. 14a), 93 days (Fig. 14b) and 140 days (Fig. 14c). It can be seen from the contour plot variable D er that the erosion indeed starts near the location of the gap, and subsequently extends along the bottom part of the pipe in depth direction. At 140 days the erosion below the pipe has become quite severe, and has induced a void with a depth of approximately half of the pipe diameter (red color) within which the erosion degradation parameter is maximal, D er ¼ D max er ¼ 0:93. Hence, the erosion void is of the suffosion type, with its volume for 95% filled by ground water, and only for 5% by remaining soil particles. It can be further observed that almost the complete saturated sand layer at this stage has undergone some degree of erosion, in agreement with an erosion degradation parameter larger than zero and below the maximum value, 0\D er \0:93. It is further interesting to notice that the erosion profile develops much faster than in the first case study on soil piping, see Fig. 10, which clearly demonstrates that the characteristic time scale of an erosion process very much depends on the type of problem and the corresponding hygro-mechanical conditions.
The groundwater flow patterns after 47 days, 93 days and 140 days are illustrated in Figs. 15a, b and c, respectively, and show to be in agreement with the erosion profiles in the area in which the flux threshold value of q fm;cr ¼ 10 À6 m/s is exceeded. Within the erosion void the permeability of the remaining ''soil structure'' is maximal, and, in accordance with Eq. (20), relates to a value of k ¼ 2:0 Â 10 À2 m/s. Note that this value is 200 times larger than the initial permeability k 0 ¼ 10 À4 m/s of the sand material, see Table 3.
The deviatoric plastic deformations generated under the dead weight loading and the subsequent erosion process are illustrated in Fig. 16 for a time instant of 140 days. In comparison with the deviatoric plastic deformations generated under only the dead weight loading, see Fig. 9b, the deviatoric plastic deformation in the soil material directly left and right of the sewer pipe has substantially increased by more than a factor of 10 to a value of c p ¼ 5:4%. Furthermore, at the bottom-left and bottom-right of the pipe the deviatoric plastic strain has reached a maximum value of c p ¼ 14:4%: Obviously, this strong, local increase in deviatoric plastic deformation is caused by the stress Fig. 15 Time evolution of the porewater velocity field under a strong groundwater flow near a defect sewer pipe (second case study), as characterized by the spatial evolution of the norm of the flux jjq fm jj (in m/s) after a 47 days, b 93 days and c 140 days; the black color indicates the area within which the value of jjq fm jj is lower than the critical threshold value q fm;cr for erosion Fig. 17 Time evolution of the vertical displacement w (in mm) of the ground surface -measured above the centre of the sewer pipe -under a strong groundwater flow near a defect sewer pipe (second case study) for i) an elasto-plastic sand material and ii) an elastic sand material Fig. 16 Spatial variation of the deviatoric plastic strain invariant c p after the application of the dead weight loading and an erosion process of 140 days under a strong groundwater flow near a defect sewer pipe (second case study) redistribution that balances the local loss of soil resistance in the suffosion erosion void below the sewer pipe; the specific local shear failure zones in the soil next to the sewer pipe are indicated in Fig. 16 by the grey hatched areas. Figure 17 depicts the vertical surface displacement as a function of time for both an elasto-plastic and an elastic sand material. The surface displacement initially grows moderately, but after approximately 100 days starts to increase rapidly due to the formation of the erosion void below the sewer pipe. The surface displacement after 140 days equals 5.8 mm and 4.4 mm in the case of, respectively, the elasto-plastic soil model and the elastic soil model. The relative difference between these displacements is substantial, i.e., 32%, which indicates the importance of accurately modelling the constitutive behaviour of the soil in erosion simulations by means of an elasto-plastic model. Also, it is interesting to notice that the ground surface deflections in Fig. 17 are much larger than those generated under soil piping erosion, see Fig. 13.
Additional simulations not presented here have shown that a decrease of the gap width in the sewer pipe by a factor of 40 to a value of 0.5 mm leads to a similar surface deflection evolution as illustrated in Fig. 17. Essentially, the decrease in gap width reduces the flow area at the bottom of the sewer pipe, which is accompanied by an increase in the local flow velocity, such that, as a net effect, the erosion void develops in a more or less comparable fashion. Nevertheless, the value of the surface deflection at 140 days is about 11% smaller, namely 5.2 mm and 3.9 mm in the case of an elasto-plastic and an elastic sand material, respectively. Note hereby that the assumption of a small gap width of 0.5 mm ignores the fact that the larger particles in the sand material can not flow into the pipe, so that the real surface deflections for this case will be smaller.
Finally, it needs to be emphasized that the remaining stability of the suffosion type of erosion void illustrated in Fig. 14c is limited, as its bearing capacity is (virtually) generated by the fluidic mixture inside; hence, as soon as the groundwater table lowers and the fluidic mixture leaves the erosion void, the sewer pipe drops down into the erosion void, whereby the overlying soil structure is likely to collapse into a sinkhole. Accordingly, the abrupt increase in surface settlement registered after 100 days can be considered as a critical warning for catastrophic failure. Since the subsequent surface displacements are in the order of several millimeters, in practice it should be possible to detect these changes with satellite radar interferometry [15,40], which may help to prevent the eventual collapse into a sinkhole.

Concluding remarks
A coupled hydro-mechanical model has been presented that describes the process of subsurface soil erosion. The governing equations for subsurface erosion have been developed by considering a saturated porous medium composed of soil particles and a fluidic mixture, and letting the particle volume fraction decrease under the development of erosion. The kinetic law proposed for the erosion process has a similar form as the type of threshold law used in erosion models at fluid-soil interfaces in channel and pipe flows. The degradation of the particle structure under erosion reduces the effective elastic stiffness of the porous medium and increases its permeability. The stiffness reduction by erosion has been incorporated in an elastoplastic soil model that describes the development of plastic deformations under frictional sliding and granular compaction.
The numerical analyses on soil piping and erosion void formation illustrate that the numerical model realistically predicts the size, location and characteristic time scale of the generated erosion profiles and the deformations of the surrounding soil structure. Hence, the erosion profiles computed by the model may be used as input for a detailed analysis of the local, residual bearing capacity and stress redistribution of buried concrete pipe systems. Additionally, the modelling results may support the early detection of in situ subsurface erosion phenomena from ground surface deformations recorded with satellite radar interferometry. The analyses further show that the characteristic times scale of the erosion process and the generated ground surface deflections strongly depend on the geometry and features of the problem, as reflected by the specific hydromechanical conditions. In addition, for a sand material the characteristic time scale for erosion generally is several orders of magnitude larger than that for soil consolidation.
In order to rigorously translate the computational results towards practical recommendations and guidelines, the material parameters of the erosion model need to be calibrated in an accurate fashion by means of a systematic experimental study, such as presented in [70]. A way to do this is to measure quantities such as the fluid flux and the amount of eroded mass in basic erosion experiments, and use these measures in an inverse way for calibrating the constitutive parameters of the erosion model. Further, the sensitivity of subsurface erosion phenomena to ground water conditions, soil type and the geometrical properties of a structural configuration are topics for future studies. Finally, the practical case studies discussed in this communication refer to two-dimensional (plane-strain) FEM models, which obviously ignore effects caused by hydromechanical fluctuations in the out-of-plane direction of the soil structure. It is worthwhile exploring these effects in detail by means of three-dimensional FEM models.

A Executing a hydro-mechanical analyis with a thermo-mechanical FEM module
The coupled hydro-mechanical analyses presented in this paper were carried out by applying the thermal-mechanical module in ABAQUS, and making use of the analogy between the partial differential equations for pore fluid transport and heat conduction. By adequately relating the model parameters of the two physical processes, the temperature h and heat flux q th following from the thermalmechanical analysis may be interpreted as the pore pressure p fm and the flux q fm of the fluidic mixture, respectively.
For establishing the relations between the process parameters of the hydro-mechanical and thermo-mechanical models, consider first the partial differential equation for the pore pressure field in a poro-elasto-plastic medium with erosion, as given by Eq. (60): Here, k is the permeability of the porous medium in accordance with Eq. (20), q s and q fm are the densities of the solid and fluidic mixture phases, _ p is the time rate of change of the total hydrostatic stress, _ D er is the rate of the erosion degradation parameter presented in Eq. (28), D er is the erosion degradation parameter given by Eq. (18), e e vol is the elastic volumetric deformation, _ e p vol is the plastic volumetric strain rate, _ m er is the rate of eroded mass defined by Eq. (16), and B is the elastic bulk modulus, which reads in which E is the Young's modulus and m is the Poisson's ratio. In addition, the partial differential equation for the process of heat conduction is [52] qc _ h À $ Á k th $h À Á ¼ r ; with q the material density, c the specific heat capacity per unit mass, k th the thermal conductivity, and r the body heat source per unit volume.
In order to further elaborate these relations, the constitutive equation for the transport of the fluidic mixture through the pores given by Eq. (19), i.e., Darcy's law: is compared to the constitutive equation for heat conduction, i.e., Fourier's law [52]: whereby the minus sign in the right-hand side of Eq. (78) indicates that the thermal flux is positive in the spatial direction of a temperature decrease. From Eqs. (77) and (78), it becomes clear that the heat flux q th can be interpreted as the flux q fm of the fluidic mixture if In addition to the above relations, the coefficient of thermal expansion a used in the thermo-mechanical analysis needs to be expressed in terms of constitutive parameters of the hydro-mechanical model. For this purpose, consider the constitutive equation for a thermo-mechanical material with erosion, as expressed in terms of the hydrostatic stress p by: whereby the final result in the right-hand side is obtained by using the strain decomposition corresponding to a thermo-elasto-plastic material e ¼ e e þ e p þ e th ; ð83Þ and the fact that the thermal volumetric strain can be written as with a the coefficient of thermal expansion. Note that in the above expression the reference temperature associated with the thermal strain is taken as zero for simplicity. Additionally, in the hydro-mechanical model the total hydrostatic stress p is decomposed in accordance with Terzaghi's stress decomposition, Eq. (17), as Equating Eq. (85) to Eq. (83) then leads to the expressions p 0 ¼ ð1 À D er ÞB e vol À e p vol À Á ; It can be confirmed that Eq. (86) 1 is in agreement with Eq. (55), and that in Eq. (86) 2 the temperature h corresponds to the pressure p fm in the fluidic mixture if whereby the expression for the bulk modulus B, Eq. (74), has been substituted to obtain the final result.
In conclusion, the temperature h and heat flux q th computed from a thermal-mechanical analysis may be interpreted as the pore pressure p fm and the flux q fm of the fluidic mixture in the hydro-mechanical erosion model if the thermal model parameters follow the expressions given by Eqs. (79), (80), (81) and (87).