Numerical analysis of failure mechanism observed in backfills supported by masonry walls

Abstract Masonry retaining walls are designed to resist lateral forces. Their stability is essentially warranted by the correct determination of the failure surface geometry. Accordingly, this study intended to investigate the influence of wall and backfill properties that control failure surface geometry of cohesionless backfills. For this purpose, the discrete element method (DEM) is utilized, and a series of parametric studies were conducted. As the wall-joint parameters reflect the mortar quality of the blocks that constitute the masonry wall, three binder types from weak to strong were defined. Additionally, loose to dense backfill soil conditions and wall-backfill interface properties were also investigated. The results indicate that in the case of a thin rigid wall, the failure surface of dense backfill is identical with the classical earth pressure theory. However, for the masonry walls with a higher foundation width, the failure surfaces are much deeper and wider; particularly on the active side compared to the classical earth pressure theories. In addition to that the deformation mechanism and the associated failure surfaces are greatly influenced from the mortar quality which results with either a deep-seated or sliding type of failure.


INTRODUCTION
Despite their widespread use, retaining walls have increased professional and academic interest in recent decades. More specifically, for masonry retaining walls the experience gained from the field observations and reassessment of the existing walls improved the perspective on the design. Masonry retaining walls are bound with mortar or it can be built with dry joints and their design must satisfy both ultimate and serviceability limit states. As it is known, external stability depends on the engineering properties of the soil, whereas internal stability is controlled by the structural performance. Accordingly, the design of retaining walls requires the consideration of soil-structure interaction. Frictional resistance and overlapping forces along the blocks, which are primarily influenced by their geometry, control the behavior and capacity of masonry retaining walls. As a result, the response of masonry retaining walls is significantly influenced by the materials used and the quality of construction, which makes capacity estimations difficult. In terms of the soil behavior, the active failure mechanism is influenced by the mode of wall movement, cohesion and internal friction angle, soil particle characteristics, interface friction and adhesion . However, determining the specific impact of each parameter in experimental investigations is difficult. As a result, numerical methods can help researchers gain a better knowledge and insight of the micro Since the external stability of a masonry wall depends on the engineering properties of the soil, the estimation of the lateral pressures acting on the wall has a critical influence on the design. The classical method suggested by Coulomb (1776) and Rankine (1857) are commonly used by practitioners for the determination of lateral earth thrust. However, these theories are based on the limit state of soils and assume that the failure plane formed in the backfill at the ultimate state is a straight line. Failure surface geometry behind the wall was examined experimentally (Leśniewska & Mróz 2000, Tejchman 2004, Leśniewska et al. 2012, Altunbas et al. 2017) by many researchers. Most of these researches are focused on the mechanical properties of the backfill, backfill density (Nuebel 2002), grain shape and diameter (Leśniewska et al. 2012, Soltanbeigi et al. 2018) dilatancy angle (Altunbas et al. 2017, Soltanbeigi et al. 2020). According to the results of these studies, the geometry of the potential failure surface in the non-limit states plays an important role for the soil response and it is not linear.
The internal stability or the structural behavior of masonry retaining walls is an active research area, especially in earthquake prone regions. Zhang et al. (2004) carried out a study to investigate the stability of drystone masonry retaining walls through equivalent continuum model and a joint element model for the wall. It is stated that more realistic behavior of a fractured medium can be simulated if joint elements are used. Alejano et al. (2012) stated that the stability of the retaining wall can be increased by forcing the value of the slope of the base part of the wall tending to overturn to be smaller. Pulatsu et al. (2021) performed static analyses to assess the capacity of the masonry retaining wall under uniform surcharge pressure using mixed discrete-continuum approach.
The Discrete Element Method (DEM), which was first introduced by Cundall (1971), has the advantage of taking into consideration; structural response, soil-wall interface, and the ability to simulate soil behavior. These advantages are critical for the internal and external stability calculations which makes 3DEC a valuable tool for the numerical modelling of masonry walls (Saygili & Lemos 2021). Additionally, 3DEC shows the movements of the jointed masses that represent the discontinuous medium as an assemblage of discrete blocks (Saygili 2019, Saygili & Lemos 2020, Saygili & Polat 2021. 3DEC uses solid physical and mathematical theory to explicitly address the geometry of topography and joints, block contacts, and block displacements. The solution is based on the integration of the equations of motion of the blocks, using an explicit time-stepping algorithm. Since the stability of the masonry type retaining wall depends on mortar properties and wall weight, the friction and cohesion between the blocks were changed to understand their effects on soil behavior.
It can be noticed that both structural and soil behavior have been considered in very limited studies (Villemus et al. 2007, Colas et al. 2010, Quezada et al. 2016, Pulatsu et al. 2020. In this work, in order to consider the associated soil and structure behavior together, DEM is used to investigate the response of a cohesionless backfill yielding to failure for different soil and structure properties. To understand the progressive response of the backfill and wall, a parametric analysis was performed. A series of analysis was carried out on the varying mechanical properties of the retaining wall and cohesionless backfill. The masonry retaining wall was represented as rigid blocks interacting with each other. 3DEC allows large displacements and rotations of blocks, thus being ADLEN ALTUNBAS et al.

NUMERICAL ANALYSIS OF FAILURE MECHANISM
able to simulate structural collapse of the wall. These blocks were lumped at the joints and joint properties that governed the structural behavior were changed with assigned contact properties. For the masonry retaining wall three different mechanical properties were assigned. To observe the effect of structural response on the failure mechanism, different joint normal and shear stiffness properties were assigned. For each case of the retaining wall model, backfill properties were changed by considering loose, medium-dense and dense sand conditions. The results clearly showed that the mechanical properties of the masonry wall greatly influence the geometry of the failure surfaces particularly on the active side as well as the deformation mode of the entire wall. The findings of this study may be useful for future research since they provide theoretical background for stability analysis and collapse prevention for similar masonry retaining walls.

Basic aspects of 3DEC
3DEC simulates the behavior of discontinuous medium i.e. jointed rock layers or masonry structures that are under static or dynamic loading. The numerical formulation is based on the discrete element method (DEM) for discontinuum modeling, and the discontinuous materials are represented as discrete blocks. The discontinuities between each block are treated as boundaries where displacements and rotations of the blocks are allowed. Each block may behave as either rigid or deformable material based on the constitutive model used to define the behavior (Itasca 2020).
The first stage in 3DEC is to establish the system's geometry, which allows the rigid blocks to dislocate freely when they come into contact. In 3DEC, computations are repeated using Newton's second law, which is applied to the joints, and the force-displacement law, which is applied to the contacts. The force-displacement law connects the joint contact forces to the displacements that arise. At each time step, the total normal and shear forces on the joints are computed by accounting for both contact and body forces acting on the block. After determining the contact forces at each joint, Newton's second law of motion is applied. It provides the acceleration, velocity, and displacements of each element between two successive time-steps. For integrating the block's equations of motion, 3DEC depends on an explicit time-stepping technique. The three translational equations of motion for a rigid block's center of mass are as follows: where u i is the displacement vector of the block center; m, the block mass; and the mass-proportional viscous damping parameter, ∝, which reproduces the energy losses in the system beyond frictional dissipation. The elastic forces are included in the force vector f i , which is the total of the applied forces, including self-weight, and the contact forces, which are a function of the relative block displacements (Saygili & Lemos 2020). Euler's equations control the three rotating degrees of freedom, which are represented in the block's principal axes of inertia as: ADLEN ALTUNBAS et al.

NUMERICAL ANALYSIS OF FAILURE MECHANISM
where i is the rotational velocity vector; I i , the principal mass moments of inertia. The total of the moments generated by the contact forces and the applied forces is the moment m i . The deformable block logic is analogous to the formulation for two-dimensional blocks, as given by Lemos (1987). All material models for deformable blocks in 3DEC assume an isotropic material behavior in the elastic range described by two elastic constants (bulk modulus, K, and shear modulus, G). Because bulk and shear moduli are thought to correlate to more basic elements of material behavior than Young's modulus and Poisson's ratio, the elastic constants, K and G, are employed in 3DEC instead of Young's modulus, E, and Poisson's ratio, . 3DEC run's solution time is determined by the number gridpoints in deformable blocks, as well as the number of contacts in a model. If the model has a small number of contacts, the time is proportional to N 3/2 , where N is the number of gridpoints in deformable blocks. Fundamentally, deformable blocks are discretized into finite-difference tetrahedral elements. Gridpoints are the vertices of the tetrahedral elements, and the equations of motion for each gridpoint are derived using (Itasca 2020): where is the mass-per-unit volume of the medium, b is the body force per unit mass, and dv/dt is the material derivative of the velocity. The local nonviscous damping used in the analyses. The damped equations of motion is: where symbol F <l> i is used to represent the sum of the contributions at global node of all tetrahedra meeting at that node, F <l> i is the damping force. n n is the total number of nodes involved in the medium representation and M <l> is the nodal mass. Local damping was used to equilibrate static simulations.

Factor of Safety Analyses -Strength Reduction Method (SRM)
In this paper, the strength reduction method is used to calculate the location of failure surfaces mobilized within the backfill behind the masonry wall. In geotechnical design, SRM has been used for decades to determine the failure surfaces and their factor of safety (Zienkiewicz et al. 1975, Griffiths 1999, Cai & Ugai 2000, Cheng et al. 2007, Wei & Cheng 2009, Kelesoglu 2016. The factor of safety (FoS) calculated by SRM is the number by which the shear strength parameters cohesion (c ′ ) and shear strength angle ( ′ ) must be factored to bring the soil mass into a state of limit equilibrium or failure. The mobilized shear strength parameters (c ′ and m ′ ) in a SRM type of failure analysis are reduced monotonically by dividing the shear strength parameters by the FoS values in small increments which are defined as: The common method to implement the SRM is to replace the mobilized c m ′ and m ′ values in place of c ′ and ′ values using an elasto−plastic soil model in a finite element or finite difference type of An Acad Bras Cienc (2023) 95(2) e20220420 4 | 18 numerical analysis. At the first step the FoS value is taken as unity, and it increases incrementally until the global failure is achieved. For 3DEC, the procedure is slightly different from the conventional finite element procedure where a characteristic response time (Nr) which is a representative number of steps that defines the response time of the system is defined first. 3DEC defines Nr by setting the strength parameters to a large value and making a large change to the internal stresses and afterwards the number of steps to return the system into equilibrium is calculated. 3DEC manual (Itasca 2020) recommends setting a maximum limit of 50,000 for Nr, where this default value is used through the analyses carried out in this paper. By default, the unbalanced force ratio is also taken as 10 −3 .

Numerical model
A series of numerical analyses were performed by varying: (i) the mechanical properties of the backfill (ii) soil wall interface friction and (iii) the stiffness of masonry retaining wall. This parametric analysis was conducted for a masonry retaining wall with a height of 8.0 m and foundation width of 5 m. Cross-section of the modelled wall is given in Fig. 1. To eliminate the effect of lateral boundaries, the horizontal dimensions of the model were chosen as 70 m. Structural components of masonry retaining wall were simulated via rigid discrete blocks, whereas the cohesionless backfill is modelled by a deformable continuum and the properties of backfill is given in Table I. Normal and shear contact stiffnesses at the interaction surface between masonry blocks are given in Table II. All system deformation is lumped at the joints in rigid block models. The joint normal and shear stiffness characterize the wall deformability. These are determined by calculating average joint spacing in each direction and multiplying it by the deformability of the mortar and block material. The normal and shear stresses in the joints are related to the relative displacements between the blocks in the normal and shear directions by joint stiffnesses. The friction between soil and masonry wall is assumed as 2/3 of the soil internal friction angle and additionally no cohesion and tensile strength assigned to this surface (Pulatsu et al. 2020). The backfill and wall parameters used in the analyses are shown in Tables I, II and III. For all these cases, the properties of the bearing stratum were not changed and modelled as a dense cohesionless foundation material. Additionally, this foundation part was replicated via a deformable block, which is represented in a fixed condition. The cohesionless backfill was modelled for three different combinations: loose, medium dense and dense conditions of sand. An elastic-perfectly plastic Mohr-Coulomb failure model with associated flow rule was defined for the backfill behavior. All the analyses were carried out in plane strain conditions. The horizontal degrees of freedom at the boundary surface of the modelled backfill are restricted whereas they are left free in the vertical direction.
Each analysis consists of two separate stages. In the first stage of the numerical analysis, the equilibrium of the model is achieved under its self-weight. At the end of this phase, FoS analysis was conducted by using the SRM. From the calculated results of the SRM, the failure surface behind and under the masonry wall is estimated using the maximum shear strain increment and velocity maps. Through the analyses of these maps, the failure surface of each case was determined. For this purpose, a coordinate system is established to be able to quantify failure surface geometry. This way, it becomes possible to observe the influence of considered parameters that may affect the soil structure interaction.

Determination of Failure Surface Geometry
The calculated failure surfaces are obtained by using a simple digitization technique. The origin of the coordinate system is defined in Fig. 2 where the lateral distance of the failure surface crest on the active side is defined as Xtop and the lateral distance of the failure surface toe on the passive side is defined as Xtoe. The failure surface is defined as the hypothetic line/interface where the shear strain values tend to increase rapidly in a very short distance as seen in Fig. 2. The shear strain increment distribution of each analysis was digitized by using this methodology. By normalizing the horizontal and vertical dimensions with the wall height (Hw), a new dimensionless coordinate system is generated. Thus, the geometries of the failure surfaces can be identified from the numerical analysis with different wall joint and backfill properties.

Mobilization of Failure Surfaces (3DEC Validation)
In order to ensure that the numerical results are realistic and reliable, it is of great importance to validate the model in the first place. To this aim, the results of 3DEC analysis are compared with those of similar experimental studies available in the literature. The details of the experimental models against which the results of this study are compared can be found in Altunbas et al. (2017).
During the tests carried out by Altunbas et al. (2017), the rigid retaining wall with dense backfill was displaced away from the backfill in the lateral direction to simulate the active failure state. In order to realistically represent these experiments with 3DEC analyses, a thin rigid wall is defined. The elasticity modulus of the wall is 40 times higher than the adjacent backfill. In the 3DEC model, the rigid wall was incrementally moved away from the backfill in a translational mode until reaching failure ( Fig.  3a-b). The outcome of the experimental comparison is presented in Fig. 3c

Influence of Wall-Joint Parameter on Failure Surface
A typical masonry wall is generally formed of bricks or stones with mortar or dry joints both in vertical and horizontal directions. Realistic analyses of masonry walls with such discrete materials, necessitate the modeling of the blocks and the interface joints between these blocks. However, in practice, when a masonry type of wall is designed, the entire body of the wall is assumed to be a single volumetric mass, where the stress-strain response of the wall is defined by using the elastic fully plastic Mohr-Coulomb model. In such cases, the binding effect of the mortar between the stones is generally represented by the cohesion intercept of the MC model. These two assumptions may lead to amplifying the strength of the wall which results in erroneously calculated higher stability and lower displacements than they should be. However, the joints of a masonry wall should be considered as a parameter that may affect the deformation mechanism mobilized through the wall itself and the backfill material. Thus, a realistic analysis of a masonry wall could be achieved only if the discrete blocks and the interface between these blocks are modelled where the interaction between the wall and the backfill is considered as well. The strength and stiffness characteristics of the wall-joint parameters gain importance when the stability of the existing structures, particularly historical ones, are investigated.
In this study, the effect of the wall joint parameters on the stability of a masonry wall is considered by changing the cohesion and tension properties of the mortar using interface joints. The discrete element model is given in Fig. 4 as the location and the geometry of the wall is shown as grey bricks. Three different wall joint parameters were defined by considering the cohesion and tension properties of the wall, namely J1, J2 and J3 which represent strong, medium and weak mortar conditions that can be observed in the field. J1 type of wall joint parameters match up with relatively strong mechanical parameters (coh=750 kPa, ten=375 kPa) and J2 type of wall with joint parameters of cohesion is 500 kPa and tension is 250 kPa whereas J3 corresponds to a wall with weak joint properties i.e. coh=250 kPa, ten=125 kPa.
The mechanical properties of the backfill behind the wall have an important effect on the wall behavior as well as the wall joint parameters. Thus, the mechanical properties of the backfill material are taken into account for dense, medium dense, and loose sand conditions. The influence of both wall joint parameters and backfill densities on the deformation mechanism and the mobilized failure surfaces were investigated through the discrete element analyses. The max shear strain increment results obtained for dense, medium dense, and loose sands associated with different wall joint parameters are given in Fig. 4. Based on the results of 3DEC analyses, general shapes of the failure surfaces considering the densities and wall joint parameters were also evaluated and plotted as shown in Fig. 5. For the three joint cases of the loose backfill condition, different wall joint models yield the same shear strain distributions which correspond to the same shear failure surfaces for all three models. In other words, it is clear from Fig. 4 that there is no considerable influence of wall joints on the failure mechanism in the loose backfill condition. These analyses were also carried out for medium dense and dense backfill materials as well. For the medium dense backfill case, the distribution of the shear strain increment for J1 and J2 is very close with small changes for J3 distribution which has a minor differentiation of the horizontal distance at the top surface of the backfill. However, considering the dense state for all three joint cases, the distribution of the shear strain diminishes from very high strain values to lower ones as the wall joint parameter decreases. In addition to that, as shown in Fig. 5, the geometry of the failure surface is getting smaller/shallower with the decrease of the wall joint parameters. As a result, the shear strain distribution observed around the retaining wall is greatly influenced by the properties of the wall and backfill. Evaluating the shear strain maps and so defined failure surfaces, it is observed that as the backfill conditions get denser the influence of the wall joint parameters becomes more dominant. For the dense backfill conditions as the wall joint parameters get stronger (J1), for the wall to reach the failure, backfill soil is required to exhibit a deep-seated failure associated with lateral thrust and lateral deformation behind the wall. However, as the joint parameters are reduced (J3), the wall behaves less intact or rigid and more deformable thus the failure of the system occurs more rapidly compared to the strong wall case (J1). As the backfill conditions get looser the failure mechanism is also altered from a deep-seated failure surface to a sliding mode since the frictional resistance developed between the base of the wall and the soil is reduced. In this case, the entire wall tends to slide instead of tilting/overturning; the influence of the wall joint parameters is only marginal.
To understand the mobilization of the failure surfaces in more detail, specifically for the dense backfill condition, additional wall joint parameters were defined to represent a very weak mortar condition, namely J4, with mechanical parameters of cohesion 125kPa and tension 62.5 kPa. The wall joint parameters of J4 case were defined to represent the mechanical properties of the old or even historical walls that has very weak mortar between the stones that constitute the wall itself.
For the dense backfill conditions, the results of maximum shear strain increments, displacement magnitudes, and displacement vectors considering different wall joint parameters, are illustrated in Fig. 6. It is evident from Fig. 6 that the failure surfaces are getting less distinguishable for the wall having J3 and J4 properties, on the contrary J1 and J2 walls have clear failure surfaces. Displacement vectors also show that as the wall joint parameters are reduced from J1 to J4 to define the reduction of the mortar strength, the amount of displacement to bring the wall into failure is also reduced. For stronger mortar properties J1, the wall behaves more rigid, and higher displacements associated with deep-seated failure surfaces are observed within the dense backfill material. As the mortar properties get weaker, the accumulated shear strains are more dominant within the wall body and a structural type of shear failure is mobilized with fewer displacements. This shows the influence of wall joint parameters on the localization of shear strains in granular materials. Based on this observation, it is possible to propose that the geometry of shear planes can be considered as functions of wall joint parameters for dense backfill.
Three types of masonry retaining wall geometry were created for the dense soil conditions, using J1, J2, and J3 wall joint properties. In Fig. 7 the variation of shear strain maps is given for the same scale for the sake of compatibility. For the second wall model, the soil in front of the wall was removed to eliminate passive lateral pressure on the face of the wall. As expected, it is observed that without passive stresses, the intensity of the strains decreases within the failure wedge. Additionally, compatible results with the previous analysis were obtained showing the influence of wall joint parameters which states the diminishing distribution of the shear strains from very high strain values to smaller ones, as the wall joint parameter decreases.

The geometry of the mobilized failure surface
In the design of the retaining walls, the location and the geometry of the failure surface is critical, not only because of the global stability of the wall itself but also to define possible countermeasures in the case of an unstable condition. The location and geometry of the failure surface obtained from J1, J2, J3, and J4 types of analyses were illustrated in Fig. 8 for the deep-seated failure in order to emphasize the influence of wall joint parameters for dense backfill material. Based on the origin point at the toe of the masonry wall, the boundary points of the obtained failure surface at the active (Xtop at Fig.  8b) and passive (Xtoe at Fig. 8a) sides are given in Fig. 8c. These results clearly show that as the wall stiffness increases in conjunction with the wall joint parameters (i.e. cohesion of the wall) the failure surfaces tend to expand in the lateral direction both through active and passive sides. It can be seen from this graph that there is almost a linear relationship between the location of the boundary point of the failure surface and the wall stiffness.
From the analyses of failure surfaces, it is possible to propose that the initial density of backfill influences the shapes of failure surfaces as it is expected (Fig. 9), similar to the observations of Altunbas et al. (2017) and Gezgin et al. (2021). Based on the results, for denser backfill conditions the failure surface gets deeper behind the retaining wall resulting in deeper and wider passive failure wedges. On the contrary, as the backfill gets looser, the wall tends to slide along the base soil rather than a deep-seated failure.

Influence of wall-backfill interface friction angle
To discuss the effect of wall-backfill interface friction angle on the failure mechanism, a series of numerical analyses were conducted with different wall roughness. To model a rough wall surface, a  An Acad Bras Cienc (2023) 95(2) e20220420 14 | 18 higher interface friction angle was assigned as the interface parameter. The results of the analysis are shown in Fig. 10. According to the analyses, it is clearly seen that the changes in the geometries of the failure planes as a result of the variation in the wall-backfill interface friction angles are insignificant for this deep-seated failure surface. Therefore, the influence of wall-backfill friction can be ignored at an active failure state as suggested by Craig (2004) and Altunbas et al. (2017). Based on the numerical findings of this paper, similar comments can be postulated not only for the active part but also for the passive part of the deep-seated failure surface.

CONCLUSIONS
The aim of this study is to provide a better understanding of the soil structure interaction of the masonry wall by considering both the wall joint parameter and the initial density of the backfill. For this purpose, DEM models are constituted, and the evolution of the failure wedges for each backfill and wall model were observed. The results were analyzed and discussed for varying conditions to better understand the interactions between the properties of both masonry wall and the backfill. The main findings are as follows: The numerical analyses carried out with 3DEC for classical rigid wall and dense backfill case yields almost identical results with the experimental findings thus it is reasonable to postulate that the DEM method is a valid tool to analyze such engineering problems.
Evaluating the shear strain maps and so defined failure surfaces, it is observed that higher backfill densities is significantly influenced by the wall joint parameters. In addition to that the backfill densities affect the shear strain distribution within the backfill as well as the failure surface geometry.
It is confirmed that for loose backfill, there is no considerable influence of wall joints on the failure mechanism of the backfill, whereas initially dense backfill is significantly affected by the mechanical properties of the structure. The geometry of the failure surface is getting smaller/shallower with the decrease of the wall joint parameters that defines the mortar/binder quality between the discrete blocks.
Shear strain distribution within the failure wedge diminishes with the decrease in magnitude of joint parameters. As the binding effect of the mortar decreases the deformation mechanism of the masonry wall tends to be a sliding type of failure. On the contrary, as the wall joint parameters increase, the wall tends to behave more rigid which results with a deep-seated type of shear failure mechanism that passes well below the wall foundation.
The geometry of the failure surface on the active side of the backfill soil is dependent on the wall geometry. With dense backfill in the case a rigid thin wall, the failure surface complies well with the experimental studies available in the literature. However, for a masonry wall with a higher foundation width, the calculated failure surfaces are much deeper and wider along the active side of the backfill.
Based on the analysis wall-backfill interface friction has no influence on both the geometry and strain distribution of failure surfaces.
Based on the comparison of our analytical findings in this study, we think that this approach may be employed in the performance analysis of retaining walls of different geometries. We use this method in a static environment with very good findings. The dynamic performance assessment of retaining walls, using real ground motion, is expected to be investigated as a continuation of the current study.