Meshless Modeling of Flow Dispersion and Progressive Piping in Poroelastic Levees

Performance data on earth dams and levees continue to indicate that piping is one of the major causes of failure. Current criteria for prevention of piping in earth dams and levees have remained largely empirical. This paper aims at developing a mechanistic understanding of the conditions necessary to prevent piping and to enhance the likelihood of self-healing of cracks in levees subjected to hydrodynamic loading from astronomical and meteorological (including hurricane storm surge-induced) forces. Systematic experimental investigations are performed to evaluate erosion in finite-length cracks as a result of transient hydrodynamic loading. Here, a novel application of the localized collocation meshless method (LCMM) to the hydrodynamic and poroelastic problem is introduced to arrive at high-fidelity field solutions. Results from the LCMM numerical simulations are designed to be used as an input, along with the soil and erosion parameters obtained experimentally, to characterize progressive piping.


Introduction
Water retaining structures, such as dams, canals, and levees, are essential for managing water resources. Often, the most economical, environmentally-conscious, and sustainable construction methods require construction using locally available soils. At times, it can be possible to construct these structures economically with masonry, concrete, and alternative materials. However, it is necessary that these structures be placed on high-quality foundation soils for satisfactory performance in the event that a solid rock stratum is not available at a reasonable depth. Therefore, the safety of these hydraulic structures depends on the capacity of the soils in the body and/or the foundation of the structure to resist the seepage forces generated by the retained water. That is, the soils should be able to safely allow the flow of retained without sustaining structural damage or erosion of soils.
Often times, the seepage on the downstream sides of hydraulic structures is found to be laden with soil particles, commonly referred to as "sand boils." Excessive seepage at specific points on the downstream side is indicative of preferential flow that may be due to anomalies in construction, foundation strata, or may indicate the formation of internal cracks voids due to movement of soils under seepage forces. Formation and gradual expansion of internal voids is frequently referred to as "progressive piping." This phenomenon may proceed unnoticed as this phenomenon occurs deep within the structure. These internal voids may collapse and form surface sinkholes, and in some cases may rapidly cause breaching failure of the structure. ICOLD [1] and Foster et al. [2] conducted statistical analyses of dam failures and found that internal soil erosion is a predominant cause of failure. Furthermore, Zaslavsky and Kassiff [3], Vallejo [4], and Kakuturu and Reddi [5] have described the conditions that promote progressive piping and outlined some basic methods for their assessment.
The success of any numerical method depends on the input of appropriate material characteristics and the experimental methods that are adopted for obtaining those characteristics. In our numerical approach, we consider the levee soil as an isotropic poroelastic material, neglecting the effects of internal particle movement. Therefore, a set of three unique values of the hydraulic conductivity, elastic modulus, and Poisson's ratio are used for calculating the pore pressures and internal stresses around a typical crack in a levee. It should be pointed here that the levee soils exhibit much higher hydraulic conductivity in the two horizontal directions compared to the vertical direction, because of the rolling compaction of the soil in horizontal layers during construction. The stresses experienced by the levee soils are usually much smaller compared to the stresses experienced by foundation soils for which a lot of data is available in the literature. Most relevant stress-strain characteristics for soils can be obtained by conducting tri-axial shear tests, which; however, produce non-linear curves with initial portions being unreliable due to experimental and equipment limitations. Therefore, it is customary to take either a tangent modulus or secant modulus for the desirable range of stresses. The triaxial sear test also facilitates measurement of volume change for reliable estimation of the Poisson's ratio.

Background
The first reservoir filling is considered very crucial to the integrity of the structure, as desiccation cracks may develop between reservoir construction and reservoir filling, and the filling may result in stresses that enlarge these cracks. A similar seepage failure mechanism is reported by Zhang and Chen [6], whereby catastrophic failure occurred in the first reservoir filling of a 71 m high Chinese rockfill dam. Additionally, cracks may develop during the service life of the structure, especially during intermittent dry periods, when foundation compression or even minor earthquakes may result in differential settlement of structural components. In laboratory-controlled conditions, Lakshmikantha et al. [7] employed image analysis techniques on clay-like soils to develop predictive models for the development of desiccation cracks and understand the conditions that influence cracking.
There are several documented instances of major piping-induced failures of water retaining structures, such as that of the Teton Dam as reported by Watts et al. [8], the Chagrin River Dam in Ohio described by Evans et al. [9], among several other cases reported by ICOLD (1995) [1]. Frequent problems faced by levees along the Mississippi, Red and Sacramento rivers. In the light of the Katrina disaster in Louisiana, public concerns have been heightened over frequent levee problems along the Mississippi, Sacramento, and Red rivers. These levee problems have been documented in several instances. Mansur et al. [10] reported the presence of sand boils caused by piping along the Mississippi river. Sills et al. [11] inferred that piping was the probable reason for failure of some levees in New Orleans that were not overtopped during the Katrina hurricane. Hagerty [12,13] provides a comprehensive and visual description of conditions that result in piping and methods for identifying piping erosion amid the complex mechanisms and interactions involved in emergent seepage failures. Gattinoni and Francani [14] analyzed the evolution of piping phenomenon and developed a model for studying the slope instability triggered by piping to understand failure conditions and applied it to the 1985 Stava Valley disaster that took place in Italy.
Internal erosion and progressive piping were detected in several experimental and in situ investigations of problem-riddled levees and dams. Through exhuming operations that confirmed piping, Johansen and Eikevik [15] described the internal erosion suffered by the Jukla Dams in Norway and the rehabilitation measures implemented. Studies by Chen et al. [16], Sjödahl [17], Inazaki [18], Rucker et al. [19], and Hung et al. [20] relied upon non-invasive geophysical methods that provided qualitative evidence of progressive piping in water retaining structures. Williams [21], Holden [22], Wang and Chen [23], and Lee et al. [24] made use of chemical or temperature tracers in combination with other methods for detection of progressive piping.
Though thousands of small cracks may develop in a hydraulic structure over the course of a year, most of these will self-heal with only a very few having a likelihood of growing larger and eroding to form the progressive piping conditions that result in failure. Some of the studies for developing a mechanistic understanding of progressive piping are by Sellmeijer and Koenders [25], Ojha et al. [26,27], Kakuturu and Reddi [5], Bonelli et al. [28], El Shamy and Aydin [29], and Bonelli and Brivois [30]. Recent research investigations by Fell et al. [31], Ojkan [32], and Yi [33] have focused on the transience and evolution of progressive piping due to natural conditions. Experimental methods are described by Kakuturu and Reddi [34] for the evaluation and characterization of some soil properties that influence progressive piping and, conversely, self-healing of cracks in earth dams.
Thus, it can be summarized that the mechanics of both fluidic and solid stratum field variables such as pore water pressures, hydraulic gradients, and principal stresses in the region surrounding the tiny cracks are the prime external drivers of crack expansion, erosion and progressive piping. These factors in conjunction with soil properties, shear strength, compressibility, gradation, interlocking, particle transportability, pore water chemistry, etc., govern the fate of tiny cracks.

Experimental Investigations of Progressive Piping
The objective of experimentation is to empirically quantify the progressive piping of cracks in soil samples, namely by applying a transient force system that represents typical hurricane forces in levees [34]. The forces in this system were estimated by considering the relevant surface water hydrodynamics on the levees. Laboratory experiments were conducted on soils that typically form the body and foundation of levees. In addition to observing piping, the experiment focused on obtaining soil properties needed for the development of the mathematical model.
A design of the experimental setup is shown in Figure 1. One of the important components of this experimental setup was the computer programmable pump that created a chosen transient force system. The flow cell housed a soil sample with pre-formed cracks of given dimensions in chosen orientations. Two simple experimental features successfully implemented were: (1) Pressure transducers in and outside the flow cell to minute pressures as the initial cracks grow into large pipes, and (2) particle size analyses of effluent to observe the mass concentration and size of eroded particles. Norway and the rehabilitation measures implemented. Studies by Chen et al. [16], Sjödahl [17], Inazaki [18], Rucker et al. [19], and Hung et al. [20] relied upon non-invasive geophysical methods that provided qualitative evidence of progressive piping in water retaining structures. Williams [21], Holden [22], Wang and Chen [23], and Lee et al. [24] made use of chemical or temperature tracers in combination with other methods for detection of progressive piping.
Though thousands of small cracks may develop in a hydraulic structure over the course of a year, most of these will self-heal with only a very few having a likelihood of growing larger and eroding to form the progressive piping conditions that result in failure. Some of the studies for developing a mechanistic understanding of progressive piping are by Sellmeijer and Koenders [25], Ojha et al. [26,27], Kakuturu and Reddi [5], Bonelli et al. [28], El Shamy and Aydin [29], and Bonelli and Brivois [30]. Recent research investigations by Fell et al. [31], Ojkan [32], and Yi [33] have focused on the transience and evolution of progressive piping due to natural conditions. Experimental methods are described by Kakuturu and Reddi [34] for the evaluation and characterization of some soil properties that influence progressive piping and, conversely, self-healing of cracks in earth dams.
Thus, it can be summarized that the mechanics of both fluidic and solid stratum field variables such as pore water pressures, hydraulic gradients, and principal stresses in the region surrounding the tiny cracks are the prime external drivers of crack expansion, erosion and progressive piping. These factors in conjunction with soil properties, shear strength, compressibility, gradation, interlocking, particle transportability, pore water chemistry, etc., govern the fate of tiny cracks.

Experimental Investigations of Progressive Piping
The objective of experimentation is to empirically quantify the progressive piping of cracks in soil samples, namely by applying a transient force system that represents typical hurricane forces in levees [34]. The forces in this system were estimated by considering the relevant surface water hydrodynamics on the levees. Laboratory experiments were conducted on soils that typically form the body and foundation of levees. In addition to observing piping, the experiment focused on obtaining soil properties needed for the development of the mathematical model.
A design of the experimental setup is shown in Figure 1. One of the important components of this experimental setup was the computer programmable pump that created a chosen transient force system. The flow cell housed a soil sample with pre-formed cracks of given dimensions in chosen orientations. Two simple experimental features successfully implemented were: (1) Pressure transducers in and outside the flow cell to minute pressures as the initial cracks grow into large pipes, and (2) particle size analyses of effluent to observe the mass concentration and size of eroded particles.   Another objective of experimentation is to evaluate the soil properties required for numerical implementation of the mathematical model of progressive piping. The soil properties that are characterized are: (1) Critical hydraulic gradient that must be exceeded for particle dislodgment in the upstream part of the crack, (2) rate of particle dislodgement once the actual hydraulic gradient exceeds the critical value, (3) critical shear stress that must be exceeded for surface erosion on the walls of the crack, (4) rate of the surface erosion once the actual shear stress exceeds that critical value, (5) parameters influencing the hydraulic fracturing, and (6) particle entrapment parameter of soil downstream of the crack. Example designs of the flow cells for obtaining these parameters are shown as Figure 2. Another objective of experimentation is to evaluate the soil properties required for numerical implementation of the mathematical model of progressive piping. The soil properties that are characterized are: (1) Critical hydraulic gradient that must be exceeded for particle dislodgment in the upstream part of the crack, (2) rate of particle dislodgement once the actual hydraulic gradient exceeds the critical value, (3) critical shear stress that must be exceeded for surface erosion on the walls of the crack, (4) rate of the surface erosion once the actual shear stress exceeds that critical value, (5) parameters influencing the hydraulic fracturing, and (6) particle entrapment parameter of soil downstream of the crack. Example designs of the flow cells for obtaining these parameters are shown as Figure 2. The results of aforementioned experimentation were used in defining parameters implemented in the numerical modeling techniques employed here.

Numerical Modeling
While conventional modeling techniques, such as the finite element methods (FEM) and finite volume methods (FVM) have been developed and routinely implemented to model complex multiphysics problems, they both require significant effort in mesh generation. The term "meshless methods" refers to a class of numerical techniques that typically rely on either a global or localized interpolation on non-ordered spatial point distributions, as opposed to a connected mesh that discretizes the domain of interest. Given the independence of an ordered mesh, there has been much interest in the development of meshless techniques, which are a relative newcomer to field of numerical methods, as they offer the hope of reducing the effort devoted to model preparation [35][36][37][38][39][40][41]. The meshless approach finds its origin in classical spectral or pseudo-spectral methods [42][43][44][45][46], which are based on global orthogonal functions such as Chebyshev or Legendre polynomials. However, these methods require a regular nodal point distribution for the implementation of orthogonal polynomials. This is contrasted with meshless methods which use a nodal or point distribution that is not required to be uniform or regular without loss of accuracy in field representation, due to the fact that most meshless techniques rely on global radial-basis functions (RBF) [47][48][49][50][51][52][53][54][55][56][57]. However, global RBF interpolation poses several issues and so, localized collocation meshless methods [58][59][60] have recently been proposed to address many of these issues.
In a series of recent publications [61][62][63][64][65][66][67][68][69][70], the co-authors have developed and validated a localized collocation meshless method (LCMM) based on the radial-basis function (RBF) interpolation for modeling of coupled viscous fluid flow, heat transfer problems, and fluid-structure interaction problems. The LCMM is shown to be computationally more efficient than a comparative finite volume method (FVM) [66,67] code whilst solving the partial differential conservation field equations of fluid flow and heat transfer on a non-ordered set of points. The method has been validated and verified against benchmark problems consisting of conventional finite volume codes for several cases, including lid-driven cavity flow, flow over a cylinder, the backward-facing step problem, and others [61][62][63][64][65][66][67][68][69][70]. The results of aforementioned experimentation were used in defining parameters implemented in the numerical modeling techniques employed here.

Numerical Modeling
While conventional modeling techniques, such as the finite element methods (FEM) and finite volume methods (FVM) have been developed and routinely implemented to model complex multi-physics problems, they both require significant effort in mesh generation. The term "meshless methods" refers to a class of numerical techniques that typically rely on either a global or localized interpolation on non-ordered spatial point distributions, as opposed to a connected mesh that discretizes the domain of interest. Given the independence of an ordered mesh, there has been much interest in the development of meshless techniques, which are a relative newcomer to field of numerical methods, as they offer the hope of reducing the effort devoted to model preparation [35][36][37][38][39][40][41]. The meshless approach finds its origin in classical spectral or pseudo-spectral methods [42][43][44][45][46], which are based on global orthogonal functions such as Chebyshev or Legendre polynomials. However, these methods require a regular nodal point distribution for the implementation of orthogonal polynomials. This is contrasted with meshless methods which use a nodal or point distribution that is not required to be uniform or regular without loss of accuracy in field representation, due to the fact that most meshless techniques rely on global radial-basis functions (RBF) [47][48][49][50][51][52][53][54][55][56][57]. However, global RBF interpolation poses several issues and so, localized collocation meshless methods [58][59][60] have recently been proposed to address many of these issues.

Mathematical Model and Governing Equations
The equations that govern the flow of a fluid through a poroelastic medium are the Navier equation coupled with the Richards equations [71,72] as: Here, p is the pore pressure, φ is the porosity of the medium, κ is the permeability of the medium, and β f , ρ f , µ f are respectively the compressibility, density, and viscosity of the fluid flowing through the porous medium.
The Navier and Richards governing equations in (1) are strongly coupled (two-way) given the derivative operators on the field variables defined by both equations; the pore pressure appears on the right-hand side of the Navier equation, while the deformation field dilatation, (∇ · → u ), appears on the right-hand side of the Richards equation.
The hydraulic conductivity, which indicates the relative permeability of the porous media to flow, is defined as: The Darcy velocity of the flow through the porous medium is defined as: Therefore, the actual front velocity of the flow is given by: The front velocity of the flow, → V f , is used to track the location of the saturation time within the porous medium by adopting a volume-of-fluids (VOF) approach [73] as in two-phase flows. Thus, the location of the saturation (phreatic) line, s, at any time-level can be traced by a transport equation as: Thus, the VOF parameter, s, which indicates the saturation of the porous media by the fluid can be used to weight the value of the mechanical properties between the saturated and unsaturated values as: Where, for instance, G u , is the shear modulus of the unsaturated medium while G s is the shear modulus of the saturated medium.
Returning to the equations in (1), notice that the Navier equation is assumed to have reached steady-state while the Richards equation is assumed to be transient to account for porous media compressibility effects. These can be addressed by explicitly evolving the Richards equation in time using a first-order (Euler) time-derivative approximation, where: Fluids 2019, 4, 120 6 of 22 while updating the deformation field, → u , quasi-statically according to the Navier equation. Additionally, it is not necessary that the deformation field, → u , be updated at every time-level but rather every few time-levels, allowing one to accelerate the solution process. A simple systematic approach is adopted to determine when the deformation field needs to be updated based on quantified changes of such field.
Furthermore, it is important to note the difference in the time scale between which the pressure signal, p, propagates throughout the medium and the time scale at which the saturation front, s, moves within the medium. The difference is most easily observed in Equation (7), above, where the fluid compressibility, β f , appears in the denominator of the right-hand side. Given the relative incompressibility of the fluid, β f is usually very small, which renders the propagation of the pressure signal, p, very fast when compared to the fluid velocity, → V f , at which the saturation line, s, moves. Thus, special attention must be given when progressing Equations (5) and (7) as different time-steps, ∆t, will be implemented to advance the equations.

Localized Collocation Meshless Method Formulation
The meshless formulation begins by defining a domain, Ω, which is circumscribed by boundary, Γ, as depicted in Figure 3. A finite set of data centers, NC, are defined and is comprised of NB points on the boundary, and NI interior points within the domain. All the NC data centers will serve as collocation points for the localized expansion of the field variables in the domain. The only difference between boundary data centers and internal data centers is simply that boundary conditions are applied to the former while governing equations are applied to the latter. approach is adopted to determine when the deformation field needs to be updated based on quantified changes of such field. Furthermore, it is important to note the difference in the time scale between which the pressure signal, , propagates throughout the medium and the time scale at which the saturation front, , moves within the medium. The difference is most easily observed in Equation (7), above, where the fluid compressibility, , appears in the denominator of the right-hand side. Given the relative incompressibility of the fluid, is usually very small, which renders the propagation of the pressure signal, , very fast when compared to the fluid velocity, ⃗ , at which the saturation line, , moves. Thus, special attention must be given when progressing Equations (5) and (7) as different time-steps, Δ , will be implemented to advance the equations.

Localized Collocation Meshless Method Formulation
The meshless formulation begins by defining a domain, Ω, which is circumscribed by boundary, Γ, as depicted in Figure 3. A finite set of data centers, NC, are defined and is comprised of NB points on the boundary, and NI interior points within the domain. All the NC data centers will serve as collocation points for the localized expansion of the field variables in the domain. The only difference between boundary data centers and internal data centers is simply that boundary conditions are applied to the former while governing equations are applied to the latter. To illustrate the meshless formulation, the diffusion for a general field variable, , in a generalized coordinate system, , time, , and a general diffusion coefficient, , will be taken into consideration as the governing equation valid in the domain, Ω, as: In addition, a set of generalized boundary conditions for the variable, , on the boundary, Γ, are given by: where , , and are imposed coefficients of ( , ) that dictate the boundary condition type and constrain values.
A linear localized expansion of a group or topology of influence points, NF, around each data center is sought such that: The terms represent unknown expansion coefficients, while the terms ( ) are expansion functions defined a priori. At each time level, a different expansion will be performed; therefore, the expansion coefficients, , will vary as time progresses.
The expansion functions, ( ), are selected from the family of radial-basis functions (RBF), particularly the inverse multiquadrics RBF. The family of RBF consist of algebraic expressions To illustrate the meshless formulation, the diffusion for a general field variable, φ, in a generalized coordinate system, x, time, t, and a general diffusion coefficient, κ, will be taken into consideration as the governing equation valid in the domain, Ω, as: In addition, a set of generalized boundary conditions for the variable, φ, on the boundary, Γ, are given by:β 1 ∂φ ∂n +β 2 φ =β 3 (9) whereβ 1 ,β 2 , andβ 3 are imposed coefficients of (x, t) that dictate the boundary condition type and constrain values. A linear localized expansion of a group or topology of influence points, NF, around each data center is sought such that: The terms α j represent unknown expansion coefficients, while the terms χ j (x) are expansion functions defined a priori. At each time level, a different expansion will be performed; therefore, the expansion coefficients, α j , will vary as time progresses.
The expansion functions, χ j (x), are selected from the family of radial-basis functions (RBF), particularly the inverse multiquadrics RBF. The family of RBF consist of algebraic expressions uniquely defined in terms of the Euclidean distance, r j (x), from a general field point, x, to an expansion point, x j . Examples of these functions include the polyharmonics RBF, multiquadrics RBF, and gaussian RBF. For the examples presented in this paper, the multiquadrics RBF with (n = 1) is chosen as the expansion function and is given such that: Multiquadrics RBF: where n is a positive exponent parameter and c is a shape parameter. The effects of these parameters on the interpolation will be discussed later. In all cases, the radial distance, r j (x), is defined in the Cartesian coordinate system as: To ensure that the general field variable approximation can capture constant and linear fields exactly, the interpolation functions are modified by NP additional polynomial functions, P j (x).
To reiterate, the globally collocated RBF-based meshless methods have some drawbacks, including poor conditioning of the ensuing algebraic set of equations. This can be addressed to some extent by domain decomposition and appropriate pre-conditioning [60]. Although very promising, these techniques can also be computationally intensive. A localized expansion approach [61][62][63][64][65][66][67][68][69][70][71][72][73] reduces the burden of the more common global interpolation methods [55][56][57][58][59][60] by limiting the field variable expansion over a smaller local topology around each data center, as opposed to the entire domain, to obtain field derivatives which are then implemented in time-marching or iterative numerical schemes. This yields multiple small interpolation matrices, which are significantly more manageable than the fully-populated global interpolation matrix of standard global. Since the approach relies on expanding known values of the field variables, it is only applicable if an explicit time-marching scheme with known initial conditions or an iterative solution scheme is formulated and is inapplicable directly to steady problems. However, this is not a drawback, as time marching can always be framed as a relaxation scheme for the iterative solution of steady-state problems. Thus, the selection of an influence region or localized topology of expansion around each data center is necessary for an accurate and unencumbered solution of the expansion and is easily accomplished by a circular search around each data center. For internal topologies, the search is automated to guarantee two conditions; first, that a minimum number of points will be included in the local topology and, second, that enough points are included in all directions around the internal data centers. For the topologies existing at domain boundaries, the search must also guarantee that points on opposing boundaries or points around a re-entry corner are not included within the topology. Figure 4 shows typical collocation topologies for internal and boundary data centers describing re-entry corners and opposing boundaries. Figure 5 shows an example of the circular search to build the topology around an internal data center of a non-uniform point distribution.   The collocation of the known field variable, f, (from previous time level or iteration step) at the points within the localized topology, leads to the following in matrix-vector form: Therefore, the expansion coefficients can be determined as: where the resulting collocation matrix is given by: And the right-hand side known vector is augmented as:     The collocation of the known field variable, f, (from previous time level or iteration step) at the points within the localized topology, leads to the following in matrix-vector form: Therefore, the expansion coefficients can be determined as: where the resulting collocation matrix is given by: And the right-hand side known vector is augmented as: The collocation of the known field variable, f, (from previous time level or iteration step) at the points within the localized topology, leads to the following in matrix-vector form: Therefore, the expansion coefficients can be determined as: where the resulting collocation matrix is given by: Fluids 2019, 4, 120

of 22
And the right-hand side known vector is augmented as: To reiterate, the NF expansion functions, χ j (x), are written in terms of the inverse multiquadrics RBF defined in Equation (11). Initially, c is determined for each expansion over the different localized topologies that span the domain based on (1) the average distance between data centers within the topology and (2) the total number of points within the topology. Then, an optimization routine is employed to modify the value of the shape parameter until the resulting collocation matrix, [C], yields a condition number in the range between 10 11 and 10 12 (in double-precision), which has been documented to produce smooth derivative fields over the domain for a wide range of test functions.
It is important to mention that the collocation matrix, [C], is dependent only on the geometrical point distribution within the localized topology. Thus, the optimization of the shape parameter, c, can be done at a setup stage when the topologies are assembled and prior to the actual solution process. However, there may be instances when optimization of the shape parameter, c, may be necessary at each solution time level, for example, when adaptive refinement is performed or when sharp discontinuities in the solution field are found.
The real advantage of the localized collocation approach is found in the calculation of derivatives of the field variable at the data center, x c , of each topology. For instance, any linear differential operator, L, can be applied over the localized expansion equation as: This can be written in matrix-vector form, allowing the evaluation of the field variable derivatives at each of the data centers, x c , to be found by a simple inner product of two small vectors: {L} which can be pre-built and stored and φ which is the updated field variable distribution in the topology of the data center. The derivation of the vector, {L}, is given in the Appendix A in Equations (A1)-(A5).
To illustrate this approach, the generalized diffusion equation, Equation (8), will be marched explicitly in time at each data center, x c , using a first-order finite-differencing approximation as: where the superscript, k, denotes the time level and ∆t denotes the size of the time step. The Laplace operator on the right-hand side can be replaced by the localized RBF interpolation as: Thus, the value of the solution field, φ , at every data center, x c , at the next time level, k + 1, can be calculated from the solution field at the current time level, k, beginning with the initial field condition. This is done very efficiently, requiring only the current field solution and a simple inner product of the linear derivative operator.
The ability to estimate the field variables and their derivatives by simple inner products of pre-built vectors is precisely what makes the LCMM an attractive solution method. Notice that while the expansion coefficients need to be rebuilt at each time level, the multiquadrics RBF only need to be evaluated at the setup stage, when these vectors are being built. This provides tremendous computational advantages over global methods. First, this reduces the CPU burden of evaluating fractional powers and complicated functions at every level of a time-marching or iterative scheme. Additionally, without the need to allocate a global collocation matrix, using only very small vectors for each data center, the memory demands of this approach are minimal.
Furthermore, imposition of the generalized boundary conditions in Equation (9), at the boundary data centers, x c , can be accomplished in a similar fashion. Recalling the general boundary condition from Equation (9) and replacing the normal derivative, ∂/∂n, at the data center, x c , with a localized expansion, yields:β where the interpolation vector, {∂n}, yields the normal derivative of the field variable, φ , at the boundary topology data center, x c . Ultimately, through some algebraic manipulation of the expression in Equation (20), a simple relation can be arrived at to determine the boundary field variable at the current time level, k + 1, as: where the boundary interpolation vector, {Γ}, is composed by a combination of the normal derivative interpolation vector, {∂n}, and the boundary condition coefficients,β 1 ,β 2 , andβ 3 , and; therefore, can be pre-built at a setup stage for every boundary data center, x c . The normal derivative interpolation vector, {∂n}, may be pre-built in one of two ways. A simple approach is to express this vector as a combination of derivative vectors in all directions and multiplied by their corresponding unit normal vector components. For instance, in 3D Cartesian coordinates, this is: For more stable results, additional internal points can be generated that "shadow" each of the boundary points. The points that make up this shadow layer are defined using the normal vector at each boundary point and point into the domain. This approach is slightly more involved; however, it allows for the direct approximation of the normal derivative at each boundary data center. The shadow layer is illustrated relative to the boundary in Figure 6. Furthermore, imposition of the generalized boundary conditions in Equation (9), at the boundary data centers, , can be accomplished in a similar fashion. Recalling the general boundary condition from Equation (9) and replacing the normal derivative, ∂/ ∂ , at the data center, , with a localized expansion, yields: where the interpolation vector, { }, yields the normal derivative of the field variable, { }, at the boundary topology data center, . Ultimately, through some algebraic manipulation of the expression in Equation (20), a simple relation can be arrived at to determine the boundary field variable at the current time level, + 1, as: where the boundary interpolation vector, { }, is composed by a combination of the normal derivative interpolation vector, { }, and the boundary condition coefficients, , , and , and; therefore, can be pre-built at a setup stage for every boundary data center, .
The normal derivative interpolation vector, { }, may be pre-built in one of two ways. A simple approach is to express this vector as a combination of derivative vectors in all directions and multiplied by their corresponding unit normal vector components. For instance, in 3D Cartesian coordinates, this is: For more stable results, additional internal points can be generated that "shadow" each of the boundary points. The points that make up this shadow layer are defined using the normal vector at each boundary point and point into the domain. This approach is slightly more involved; however, it allows for the direct approximation of the normal derivative at each boundary data center. The shadow layer is illustrated relative to the boundary in Figure 6. This shadow layer approach mitigates the inaccuracies of directional derivative vectors for topologies that exist at, and are truncated by, the domain boundary. This approach is especially advantageous in corners and highly curved boundaries. With the implementation of the shadow layer, the normal derivative interpolation vector, { }, can be computed simply as: This shadow layer approach mitigates the inaccuracies of directional derivative vectors for topologies that exist at, and are truncated by, the domain boundary. This approach is especially advantageous in corners and highly curved boundaries. With the implementation of the shadow layer, the normal derivative interpolation vector, {∂n}, can be computed simply as: where, r s is the distance from the boundary data center to its corresponding shadow point in the domain. Notice that all elements of the interpolation vector except two vanish. This shadow layer approach also allows the implementation of higher-order derivatives by simply including additional shadow layers in the normal direction of the boundary data center. In addition, this approach can be implemented to estimate not only the normal derivatives but the tangential derivatives as well. That is, an interpolation in the tangential direction can be formulated such that the tangential derivative at the boundary data center, x c , is estimated as: where the tangential derivative interpolation vector, {∂s}, is found through direct finite-differencing in the tangential direction, similarly to the way it is done in the normal direction. This approach is very useful as it can be employed to estimate other directional derivatives by simple rotations. For instance, in 2D Cartesian coordinates, the x and y derivatives can be computed from the normal and tangential derivatives as: To take full advantage of the highly accurate field variable interpolation capabilities, the field variable can be interpolated with the RBF to virtual locations and passed to the data center like, for instance, to a set of locations where finite-differencing approximation of derivatives can be performed. This can be implemented in the same localized topologies previously defined for the LCMM and formulated, as it will be seen shortly, to pass the desired information to the data center via an up-winding scheme, as it will be shown in the following section.
To begin, a set of "virtual" points are distributed in the topology at the locations required by finite-differencing evaluation of the derivatives (i.e., at points n, s, e, w, indicating the north, south, east, and west directions, respectively). It is important to note that these virtual points are not additional internal data centers added to the field. Instead, they are locations at which the field will be interpolated. However, the spacing employed in the distribution of the virtual points is consistent with the average spacing of data centers within the topology. This virtual point approach is illustrated in Figure 7.
Using a finite-differencing approach, the field derivatives can be evaluated through the RBF-interpolation at the virtual locations. For example, the Laplace operator can be defined at the data center, x c , assuming uniform virtual spacing, h, in all directions, using the RBF-interpolated field value at the north, south, east, and west virtual points, x n , x s , x e , x w , respectively, such that: Of course, it is not necessary to require uniform spacing, h, of virtual points within the topology, and an equivalent form of the Laplace operator can be derived if such a case were implemented. Additionally, it is worth mentioning here that any derivative operator can similarly be defined with a finite-differencing approximation using virtual points within the topology.
instance, to a set of locations where finite-differencing approximation of derivatives can be performed. This can be implemented in the same localized topologies previously defined for the LCMM and formulated, as it will be seen shortly, to pass the desired information to the data center via an up-winding scheme, as it will be shown in the following section.
To begin, a set of "virtual" points are distributed in the topology at the locations required by finite-differencing evaluation of the derivatives (i.e., at points n, s, e, w, indicating the north, south, east, and west directions, respectively). It is important to note that these virtual points are not additional internal data centers added to the field. Instead, they are locations at which the field will be interpolated. However, the spacing employed in the distribution of the virtual points is consistent with the average spacing of data centers within the topology. This virtual point approach is illustrated in Figure 7. Using a finite-differencing approach, the field derivatives can be evaluated through the RBFinterpolation at the virtual locations. For example, the Laplace operator can be defined at the data center, , assuming uniform virtual spacing, ℎ, in all directions, using the RBF-interpolated field value at the north, south, east, and west virtual points, , , , , respectively, such that: The field variable, φ(x), can be evaluated at any of the virtual locations, for instance, x e , using the polynomial-augmented RBF interpolation expansion described in the previous sections, reducing to: Likewise, this can be written in matrix-vector form, allowing the offset evaluation by formulation of a directional interpolation vector, denoted as {I e } for the easterly direction. The derivation of this interpolation vector is given in Equations (A6)-(A10), resulting in: where the east interpolation vector, {I e }, is given by: The expression in Equation (28) can be easily extended to other Virtual points as: Therefore, the Laplace operator over the field variable, φ(x), at the data center, x c , reduces to: where: This definition still retains the capability of defining field variable derivatives at the data center, x c , through the inner product with a derivative interpolation vector, {L}, which is also pre-built and stored at the setup stage of the problem. Again, for consistency in the field interpolation within a localized topology, this can be implemented over the exact same local topology employed for the localized RBF interpolation and moving least-squares smoothing.
Note that in the case that such a virtual point coincides with an actual meshless data center within the same localized topology, the interpolation vector, {I}, at the virtual point in question, reduces to zero in all entries except for at the entry of the coincident point, where the vector evaluates to one. For instance, in the case of the interpolation vector, {I c }, where the virtual point always coincides with the actual point, x c , the interpolation vector reduces to: Similarly for the cases where any other virtual point (x e , x w , x n , x s , · · ·) coincides with an actual point in the topology. This dramatically reduces the computational time required for setup, especially in cases where the topologies contain uniformly-distributed points, as it is no longer necessary to optimize the shape parameter, c, for the topology and no matrix inversion is necessary to generate derivative interpolation vectors.
This finite-differencing enhanced LCMM approach proves advantageous in up-winding schemes for highly convective flow solutions, as will be seen in the following section, because it affords the possibility of having control over virtual point location within the topology, allowing for higher-order and offset finite-differencing approximations of the derivatives.

Up-Winding Schemes
In the cases when the field variable is not only diffused within the domain but convected by a flow field, → V, the governing equation for the field variable becomes a transport equation as: Special attention should be given to the convective derivative term, ( → V · ∇)φ, in the event that a threshold is surpassed where the convective forces dominate the diffusion forces. A grid-based Reynolds number is employed for this purpose as: where h can be assumed to be the average point spacing within the topology. The usual threshold that requires special treatment of the convective derivative is when Re h > 1, in which case an upwinded or offset derivative in the opposite direction of the flow field, → V, is necessary to achieve stability of the numerical solution scheme.
For example, in the case of 2D Cartesian coordinate system, the convective derivative has the following form: where V x and V y are the x and y components, respectively, of the flow field, → V. Concentrating on the x-direction and placing a set of virtual points in the topology as illustrated in Figure 8, three types of finite-differencing schemes for offset derivatives in the x-direction may be formulated and are given in Equations (A11)-(A13).
( ⃗ ⋅ ∇) = + (36) where and are the and components, respectively, of the flow field, ⃗ . Concentrating on the x-direction and placing a set of virtual points in the topology as illustrated in Figure 8, three types of finite-differencing schemes for offset derivatives in the x-direction may be formulated and are given in Equations (A11)-(A13). As it is well known, high-order finite-differencing approximation schemes are more accurate but tend to yield unnatural oscillations that can render the solution process unstable, while first-order finite-differencing approximation schemes are very stable at the expense of artificially diffusing the solution. Therefore, a blend between low-and high-order approximation schemes is necessary to compromise between accuracy and stability. In any case, all these operators can be pre-built and stored at the problem setup stage in the same way it is done for the standard differential operators discussed in previous sections, that is, for example: As it is well known, high-order finite-differencing approximation schemes are more accurate but tend to yield unnatural oscillations that can render the solution process unstable, while first-order finite-differencing approximation schemes are very stable at the expense of artificially diffusing the solution. Therefore, a blend between low-and high-order approximation schemes is necessary to compromise between accuracy and stability. In any case, all these operators can be pre-built and stored at the problem setup stage in the same way it is done for the standard differential operators discussed in previous sections, that is, for example: where ∂x 1 w T , ∂x 2 w T , and ∂x 3 w T are the first-, second-, and third-order x-derivative interpolation vectors, respectively, up-winded from the west to approximate the derivative at the topology data center, x c .
A limiter is introduced to the current up-winding technique to suppress potential oscillations. Compared to the cell-based finite volume method (FVM), the LCMM is a point-based method where no cell face values are reconstructed. Hence, a direct application of the flux corrected transport (FCT) [74,75], the total variation diminishing (TVD) [76][77][78][79] or normalized variable diagrams (NVDs) [80][81][82] limiter schemes is not possible. However, the FCT, TVD, and NVDs concepts of monotonicity could be used as basis to form a limiter for the LCMM. To this end, the field variable is first evaluated at the Vvirtual points (downstream point, the data center, the upstream point, and one far upstream point) to detect if a minimum or maximum occur at the data center or upstream point, meaning that the field is not monotonic on that particular direction. If neither a minimum nor a maximum is detected at these locations no limiting action is needed and major up-winding (third-order approximation of the derivative) is performed. If a minimum or a maximum is detected at any of these locations a limiting action is introduced and minor up-winding (first-order approximation of the derivative) is performed. This minor up-winding leads to stable solutions; however, these are highly dissipative. Similar up-winding techniques were reported in [83].

Dam-Breaking Problem
The progressive piping process is driven to a large extent by the pumping mechanism caused by transient hydrodynamics. To this end, a test of the LCMM numerical methodology is a transient hydrodynamics problem that simulates waves generated by tidal influence or hurricane force winds. The classical benchmark problem of the instantaneous dam-breaking in a "sloshing" tank was chosen for this purpose and solved using the LCMM volume of fluid (VOF) approach adopting a standard pressure-velocity correction scheme for the fluid flow problem [84][85][86][87]. This approach has been previously implemented and validated by the authors in highly-convective flows using the LCMM, see [70][71][72][73]. Here, the tank was initially filled with air (s = 0) everywhere except on the bottom left corner where a block of water (s = 1) was trapped by a dam, as seen in Figure 9a. The tank was 20 m wide by 10 m high and the initial water content was 8 m wide by 4 m high. The LCMM point distribution is also shown in Figure 9b.

Dam-Breaking Problem
The progressive piping process is driven to a large extent by the pumping mechanism caused by transient hydrodynamics. To this end, a test of the LCMM numerical methodology is a transient hydrodynamics problem that simulates waves generated by tidal influence or hurricane force winds. The classical benchmark problem of the instantaneous dam-breaking in a "sloshing" tank was chosen for this purpose and solved using the LCMM volume of fluid (VOF) approach adopting a standard pressure-velocity correction scheme for the fluid flow problem [84][85][86][87]. This approach has been previously implemented and validated by the authors in highly-convective flows using the LCMM, see [70][71][72][73]. Here, the tank was initially filled with air ( = 0) everywhere except on the bottom left corner where a block of water ( = 1) was trapped by a dam, as seen in Figure 9a. The tank was 20 m wide by 10 m high and the initial water content was 8 m wide by 4 m high. The LCMM point distribution is also shown in Figure 9b. The dam was instantaneously removed at = 0 and the water was allowed to slosh towards the right-hand side of the tank. Plots of the time-progression of the water sloshing is shown in Figure  10 for time values of = 3, 6, 9, 12 s. A plot of the velocity vectors is shown in Figure 11 revealing the motion of the air ( = 0) induced by the water ( = 1) sloshing. The dam was instantaneously removed at t = 0 and the water was allowed to slosh towards the right-hand side of the tank. Plots of the time-progression of the water sloshing is shown in Figure 10 for time values of t = 3, 6, 9, 12 s. A plot of the velocity vectors is shown in Figure 11 revealing the motion of the air (s = 0) induced by the water (s = 1) sloshing.

The Levee Problem
The final example deals with the problem of water flow through a compacted sand poroelastic levee. This problem combines the effects of storm surges to the levee as well as the progressive piping within the levee due to the pore pressure and stress field. The schematics of this problem are illustrated in Figure 12 Figure 13 after a steady-state solution is achieved. The time Figure 11. Velocity vector plot of the dam-breaking problem solution at t = 12 s.

The Levee Problem
The final example deals with the problem of water flow through a compacted sand poroelastic levee. This problem combines the effects of storm surges to the levee as well as the progressive piping within the levee due to the pore pressure and stress field. The schematics of this problem are illustrated in Figure 12. The trapezoidal portion of the levee is 40 m at the base and 3 m at the top with a height of 6 m. The water level on the left-hand side of the levee is 1 m from the top of the levee. A 3 m-long pipe (crack) inclined 45 • is included 8 m from the left-hand corner of the trapezoidal section of the levee. The levee properties where characterized as: hydraulic conductivity K h = 10 −7 m/s, Poisson ratio ν = 0.35, modulus of elasticity E = 120 MPa, dry density ρ = 1500 kg/m 3 , specific gravity of soil particles SG = 2.65, and porosity φ = 0.3. The water density, viscosity, and compressibility are ρ f = 1000 kg/m 3 , µ f = 10 −3 Pa · s, and β f = 4.55 · 10 −10 m 2 /N, respectively. The contour plots of hydrostatic pressure levels and water flow vectors within the levee and around the pipe are shown in Figure 13 after a steady-state solution is achieved. The time progression of the saturation front within the poroelastic levee at t = 20, 40, 60, 80, 100, 120 days is shown in Figure 14. progression of the saturation front within the poroelastic levee at = 20, 40, 60, 80, 100, 120 days is shown in Figure 14.  progression of the saturation front within the poroelastic levee at = 20, 40, 60, 80, 100, 120 days is shown in Figure 14.

Conclusions
The formation and expansion of internal voids within a water retaining structure, referred to as "progressive piping," is a predominant cause of failure in such structures. Given that crack formation and propagation is internal to the structure, and the methods required to obtain in situ measurements, progressive piping may occur undetected. It is understood that crack expansion is dependent on hydraulic and mechanical conditions in their vicinity. Therefore, fundamental consideration of these conditions would provide insight to the mechanism of progressive piping.
To this end, the localized collocation meshless method (LCMM) is applied to the numerical analysis of pore water pressures, hydraulic gradients, and principal stresses around such cracks in levees, both here and in prior work. The LCMM formulation is an adaptation to conventional

Conclusions
The formation and expansion of internal voids within a water retaining structure, referred to as "progressive piping," is a predominant cause of failure in such structures. Given that crack formation and propagation is internal to the structure, and the methods required to obtain in situ measurements, progressive piping may occur undetected. It is understood that crack expansion is dependent on hydraulic and mechanical conditions in their vicinity. Therefore, fundamental consideration of these conditions would provide insight to the mechanism of progressive piping.
To this end, the localized collocation meshless method (LCMM) is applied to the numerical analysis of pore water pressures, hydraulic gradients, and principal stresses around such cracks in levees, both here and in prior work. The LCMM formulation is an adaptation to conventional numerical techniques, leveraging localized collocation of radial basis functions to create a robust and accurate solution scheme to model poroelastic behavior through a VOF approach, while reducing the effort required to generate said model. Additionally, the LCMM scheme lends itself to the implementation of more advanced numerical formulations, such as up-winding with a limiter, to more accurately resolve the field variables in question.
The behavior of the fluid and porous media as modeled through the LCMM are used to develop a mechanistic understanding of the conditions necessary to prevent piping and to enhance the likelihood of self-healing of cracks in levees. The field results, such as principal stresses, from the enhanced LCMM approach are to be used as an input to the development of the progressive piping model. This modeling approach, used in conjunction with empirically-derived soil and erosion parameters, are used to further develop the model for progressive piping. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest. In matrix-vector form:
Substitution of the expansion coefficients, α j , found in Equation (14), leads to: where the east interpolation vector, { e }, is given by: Finite-Differencing Derivative Offset Schemes: (i) First-order: (ii) Second-order: (iii) Third-order: