Slope stability analysis of a landfill subjected to leachate recirculation and aeration considering bio-hydro coupled processes

During the operation of landfills, leachate recirculation and aeration are widely applied to accelerate the waste stabilization process. However, these strategies may induce high pore pressures in waste, thereby affecting the stability of the landfill slope. Therefore, a three-dimensional numerical analysis for landfill slope stability during leachate recirculation and aeration is performed in this study using strength reduction method. The bio-hydro coupled processes of waste are simulated by a previously reported landfill coupled model programmed on the open-source platform OpenFOAM and then incorporated into the slope stability analysis. The results show that both increasing the injection pressure for leachate recirculation and maximum anaerobic biodegradation rate will reduce the factor of safety (FS) of the landfill slope maximally by 0.32 and 0.62, respectively, due to increased pore pressures. The ignorance of both waste biodegradation and gas flow will overestimate the slope stability of an anaerobic bioreactor landfill by about 20–50%, especially when the landfilled waste is easily degradable. The FS value of an aerobic bioreactor landfill slope will show a significant reduction (maximally by 53% in this study) when the aeration pressure exceeds a critical value and this value is termed as the safe aeration pressure. This study then proposes a relationship between the safe aeration pressure and the location of the air injection screen (i.e., the horizontal distance between the top of the injection screen and the slope surface) to avoid landfill slope failure during aeration. The findings of this study can provide insights for engineers to have a better understanding of the slope stability of a bioreactor landfill and to design and control the leachate recirculation and aeration systems in landfills.


Introduction
Sanitary landfills still play a significant role in the disposal of Municipal Solid Waste (MSW) in most countries. Technologies including aeration and leachate recirculation have been applied to some landfills (Ritzkowski and Stegmann 2012;Stegmann 2019;Townsend et al. 2015) to accelerate the waste stabilization and reduce the methane emission potential. However, the addition of pressurized air and leachate would generate high pore pressures in landfilled waste and reduce its effective stress and shear strength (Byun et al. 2019;Stoltz et al. 2012). An unreasonable design of aeration and recirculation systems might even cause catastrophes such as landfill slope failure and incidental environmental pollution and casualties (Feng et al. 2021a;Koerner and Soong 2000;Lavigne et al. 2014). Thus, it is important to evaluate the influence of leachate recirculation and aeration on the landfill slope stability and provide some guidelines that balance the efficiency of waste stabilization and the safety of landfill slope.
The distribution of pore pressures and the change in landfill slope stability during leachate recirculation have been studied by many researchers. For instance, Xu et al. (2012) and Feng et al. (2018) utilized the limit equilibrium method (LEM) and strength reduction method (SRM) to evaluate the stability of a landfill slope recirculated by horizontal wells and vertical wells, respectively. Feng et al. (2018) also developed a simple design method for vertical leachate injection wells considering both slope stability and recirculation efficiency. However, the waste biodegradation and landfill gas flow that would significantly affect the redistribution of pore pressures and the reduction of waste effective stress during recirculation (Lu et al. 2019) were generally neglected in most numerical analyses for landfill slope stability (Byun et al. 2019;Feng et al. 2018;Koerner and Soong 2000;Xu et al. 2012). Compared with soil, MSW is characterized by complicated biodegradation that is sensitive to the environment such as water saturation and oxygen concentration. Many environmental disasters induced by landfills are closely related to the biodegradation of MSW, and the generation and migration of gas and leachate in landfills . Therefore, the stability analysis of a bioreactor landfill slope should consider waste biodegradation and the corresponding multi-phase flow.
For landfill aeration, several numerical models under different conditions have been developed by previous researchers (Cao et al. 2018;Feng et al. 2019Feng et al. , 2021bFytanidis and Voudrias 2014;Ma et al. 2020;Omar and Rohani 2017). For example, Fytanidis and Voudrias (2014) considered both the waste aerobic biodegradation and multi-phase flow in landfills under aeration based on the finite volume method (FVM). Cao et al. (2018) further considered the existence of both the aerobic reaction and anaerobic reaction of waste in a landfill.  presented a comprehensive overview of the mathematical formulations of biochemical, hydraulic, mechanical, and thermal processes in this coupled problem. These studies make a significant contribution to better understanding the coupled processes in aerobic landfills. However, to the best of the authors' knowledge, few studies have further focused on the impact of aeration on the landfill slope stability, which might be a concern for landfill engineers.
In this paper, a three-dimensional numerical analysis for landfill slope stability during leachate recirculation and aeration is performed using strength reduction method. Different from previous studies, the redistribution of pore pressures in waste is calculated by a previously reported two-phase flow model with fluid generation due to anaerobic/aerobic waste biodegradation. Hence, the responses of accelerated waste biodegradation and leachate-gas flow during recirculation and aeration can be incorporated into the slope stability analysis. Various recirculation and aeration scenarios are then simulated to provide some guidelines for the design of leachate/air injection systems considering the slope stability of landfills.

Methods
This section firstly introduces the bio-hydro coupled model for waste that was reported and validated by . On this basis, the stability analysis method for landfill slope and the numerical implementation of the entire process are presented in this study to investigate the effects of leachate recirculation and aeration. The following assumptions were adopted in the numerical modelling: (1) landfill gas and leachate were immiscible, and their migrations can be described by Darcy's law; (2) leachate was incompressible, and landfill gas was assumed as ideal gas; (3) the intermediate products of waste biodegradation processes were neglected; (4) and the small-strain assumption was applied to calculate waste deformation.

Bio-hydro coupled model Fluid flow equations
Gas and leachate flow in landfills can be described by the classic two-phase flow model for unsaturated porous materials as follows: where n is the porosity of waste; S α , ρ α (kg/m 3 ) and v α (m/s) are the saturation degree, density, and Darcy velocity of phase α (g for gas and l for leachate), respectively, and S g + S l = 1; and Q α (kg/m 3 /s) is the source term of phase α due to biodegradation. The Darcy velocities of gas and leachate can be calculated according to Darcy's law as: where the tensor field K α (m 2 ) is the permeability of phase α in waste; μ α (kg/m/s) and p α (Pa) are the dynamic viscosity and pore pressures of phase α, respectively; and g (m/s 2 ) is the gravitational acceleration.
The leachate and gas permeability of landfilled waste can be expressed as: where A is the anisotropic coefficient of waste; K v (m 2 ) is the vertical intrinsic permeability of waste; and k rα is the relative permeability of phase α, which can be estimated by van Genuchten-Mualem model as follows (Mualem 1976;van Genuchten 1980): where S le , S m , and S r are the effective, maximum and residual leachate saturations, respectively; and m is a dimensionless constant for the model. The relationship between S le and pore pressures (p l and p g ) can be expressed by van Genuchten model as van Genuchten 1980): where p c (Pa) and p c0 (Pa) are the capillary pressure and the entry capillary pressure of gas, respectively. The density of leachate was assumed as constant, while the density of gas mainly depends on the gas pressure, and its time derivative of density can be written as: where M g (g/mol) is the average molecular weight of gas mixture; and T (K) is the temperature.Substituting Eqs.

Oxygen transport equations
In aeration scenarios, the oxygen concentration can significantly affect the reaction modes and the degradation rate of waste (Kim et al. 2007;Omar and Rohani 2017), thereby affecting the distribution of pore pressures. To evaluate the oxygen distribution in landfills, its mass conservation in gas phase can be written as: where C O (kg/m 3 ), J O (kg/m 2 /s), and Q O (kg/m 3 /s) are the concentration, diffusive flux and source term of oxygen in gas phase, respectively. J O can be calculated using Fick's law: where D O (m 2 /s) is the diffusion coefficient of oxygen in gas; and τ is the tortuosity factor for gas diffusion considering the effect of porosity and saturation, and is defined as (Millington and Quirk 1961): The time derivative of C O can be written as: where Y O is the mass fraction of oxygen ( . Incorporating Eq. (12) into Eq. (9) yields:

Biodegradation equations
The reaction mode of waste in landfills depends on the local concentration of oxygen. Kim et al. (2007) suggested a threshold pressure of oxygen required for aerobic biodegradation: where p Og (Pa) is the oxygen partial pressure and RM is the parameter used to represent different reaction modes. To calculate the source term due to biodegradation mentioned in Eqs. (7), (8) and (13), a Monod-type biodegradation sub-model considering the effect of leachate saturation was adopted in this study (El-Fadel et al. 1996;Fytanidis and Voudrias 2014;Omar and Rohani 2017), which can be expressed as: where R A (kg/m 3 /day) and R N (kg/m 3 /day) are the growth rates of aerobic and anaerobic species, respectively; X A (kg/m 3 ) and X N (kg/m 3 ) are the concentration of aerobic and anaerobic species, respectively; k A,max (day −1 ) and k N,max (day −1 ) are the maximum biodegradation rates under aerobic and anaerobic conditions, respectively; f s is the leachate saturation correction factor; S (kg/m 3 ) is the biodegradable substrate concentration; S A (kg/m 3 ) and S N (kg/m 3 ) are the half-saturation constants of the substrates for aerobic and anaerobic species, respectively; k O (kg/ m 3 ) is the oxygen half-saturation constant; and R A,D (kg/ m 3 /day) and R N,D (kg/m 3 /day) are the decay rates of aerobic and anaerobic species, respectively. R A,D and R N,D can be expressed as: where X N,0 (kg/m 3 ) and X A,0 (kg/m 3 ) are the initial concentrations of anaerobic and aerobic species, respectively. The leachate saturation correction factor f s can be expressed as (Meima et al. 2008): The process of waste biodegradation can be described using the following chemical equations (Tchobanoglous et al. 1993): where C a H b O c is the molecular formula of waste depending on its chemical composition; and a, b and c are the constants representing the contents of carbon, hydrogen and oxygen, respectively (Feng et al. 2021b). According to the law of mass conservation and Eq. (20), the production rates of methane R M (kg/m 3 /day) and carbon dioxide R C (kg/m 3 /day), and the consumption rates of oxygen R O (kg/m 3 /day) and water R H (kg/m 3 /day) can be calculated as follows: where Y N and Y A are the biomass/substrate yield coefficients of anaerobic and aerobic conditions respectively. Therefore, the aforementioned source terms Q l , Q g , and Q O can be written as:

Slope stability analysis
The finite volume method (FVM) is usually used to deal with computational fluid dynamics problems. The introduced bio-hydro coupled model for waste has been programmed into the OpenFOAM platform based on FVM by . This section is to add a sub-program for slope stability analysis based on the strength reduction method, which can expand the application of FVM in the deformation analysis of solid phase. The pore pressures calculated by the bio-hydro coupled model will be directly imported into the slope stability analysis. The effective stress can be calculated based on the effective stress principle of unsaturated soil (Khoei and Mohammadnejad 2011): where σ′ (Pa) and σ (Pa) denote the effective stress and total stress, respectively; b is the Biot coefficient; I is the identity tensor; and p (Pa) is the average pore pressure, which can be expressed as: This study adopts the Mohr-Coulomb model to calculate the mechanical plastic strain of waste, as has been applied in several studies Lu et al. 2019: where f and g are the yield function and plastic potential function, respectively; σ 1 ′ (Pa) and σ 3 ′ (Pa) are the maximum and minimum principal stresses, respectively; and c (Pa), φ (°), and ψ (°) are the cohesion, friction angle and dilation angle of waste, respectively. The plastic strain increment can be expressed as: where Λ is the plastic factor that can be calculated using the local return mapping method proposed by Clausen et al. (2007).
Based on the small-strain theory,  derived the relationship between the plastic strain increment and deformation increment of waste as: where μ and λ are Lame constants which can be derived from Young's modulus E s and Poisson's ratio ν; and u (m) is the displacement vector.
This study adopts the strength reduction method (SRM) to analyze the stability of a landfill slope during recirculation and aeration. The SRM is suitable for linear Mohr-Coulomb failure criterion and has been widely used in the stability analysis of rock slopes (Yuan et al. 2020), soil slopes (Dawson et al. 1999;Griffiths and Lane 1999), and landfill slopes . The factor of safety is defined as the number by which the original shear strength parameters are divided to bring the slope to the point of failure, which is the same as that used in the traditional limit equilibrium method (Griffiths and Lane 1999). Therefore, the response of waste deformation can be examined by reducing the shear strength parameters in Eq. (28) following: where FS is the factor of safety. By gradually increasing the value of FS, the slope will eventually reach an unstable state at which the displacement calculated by Eqs. (26)-(30) shows a dramatic increase. The FS at this moment reaches the real FS.

Numerical model implementation
Equations (7), (8), (13), and (30) are the main governing equations for the landfill coupled model established in this study. There are four independent unknown variables: leachate pressure p l , gas pressure p g , mass fraction of oxygen in gas phase Y O , and displacement increment du. Other variables can be derived based on the above basic variables.
Based on FVM, the governing equations of the coupled model are discretized with the Gauss divergence theorem which converts the volume integrals of the governing equations over a specific volume into surface integrals. The detailed information on the discretization process can be found in  and Weller and Tabor (1998). Finally, the implicit term is converted into a matrix of unknown variables; and the explicit terms such as source terms can be calculated using the results obtained in the previous iteration or time step.
After discretizing, the governing equations in the centroid C of each control volume can be transformed and rearranged into linear algebraic equations as follows: where β C , β N , and R C are the diagonal coefficients, neighbor coefficients, and source terms, respectively. Assembling Eq. (32) for all control volumes as a whole, a system of linear equations can be obtained as: where B is a symmetric matrix; F is the unknown vector for leachate pressure, gas pressure, mass fraction of oxygen, and displacement increment; and b is the right-hand side vector. Figure 1 shows the overall procedures of the solver based on the sequential iteration method. At each iteration of time step t, the fluid flow equations, the oxygen transport equation, and the displacement equation are solved in sequence, with the related parameters and source terms being updated. The convergence criteria of the calculation are defined as: where ε F and ζ F are the absolute and relative convergence tolerances for variable F, respectively.The local return mapping method is adopted to obtain the modified plastic strain increment using the pore pressure calculated in the previous section. The detailed information about the algorithm can be found in Clausen et al. (2007) and Tang et al. (2015). The strength parameters c and φ are gradually reduced to examine the response of deformation until a dramatic increase in the displacement. Based on FVM, the code for the waste bio-hydro coupled model and landfill slope stability analysis are implemented in the C++ based open-source platform OpenFoam. The code is flexible in dimensionality and can solve problems in one dimension, two dimensions, and three dimensions (Lu et al. 2019.

Conceptual model and parameter selection
A three-dimensional conceptual model of a landfill slope is established with a cover on the top and a leachate collection and removal system at the bottom (Fig. 2a). In leachate recirculation scenarios (Fig. 2b), a horizontal trench (1 m × 1 m) for leachate addition with a continual pressure of P li (kPa) is installed 30 m above the bottom of a 50-m high and 180-m wide landfill with a slope gradient of 1:λ = 1:3. The slope geometry of leachate recirculation scenarios is the same as those in Xu et al. (2012) and Feng et al. (2018) for comparison of FS in the next section. In aeration scenarios (Fig. 2c), the landfill slope has a height of 25 m, a width of 80 m and a slope gradient of 1:λ = 1:3. A smaller landfill height is adopted because the influence radius of aeration wells is generally less than 10 m (Fytanidis and Voudrias 2014), and hence the effect of aeration on the slope stability can be presented more clearly. A vertical air injection screen with a length of 3 m, a diameter of 0.3 m, and an injection pressure of P gi (kPa) is placed 15 m vertically above the bottom of the landfill. The horizontal distance between the top of air injection screen and the slope surface is defined as d. The d ranges between 5 and 10 m in the following to investigate its effect on FS.  Table 1. The surrounding boundaries (the left, back, and front boundaries in Fig. 1a) were set as symmetrical. The top boundary can move freely during waste deformation and was regarded as impermeable to leachate without consideration of rainfall infiltration and evaporation. An atmospheric pressure was also set at the top. At the bottom surface, the leachate pressure was fixed at zero for free drainage and the displacement was constrained in every direction. Both the bottom surface and the leachate injection trenches were assumed as impermeable to gas since the surrounding waste is generally saturated by leachate. Similarly, the air injection screens were assumed as impermeable to leachate and the mass fraction of oxygen was fixed at 26.6% as in the atmosphere. The right surface was assumed as impermeable to both leachate and gas due to the existence of an engineered berm at the toe of landfill slope. The initial stress and displacement fields are obtained by solving Eqs. (28)-(30) under gravity.
Default model parameters for slope stability analysis are listed in Tables 2, 3 and 4. The biokinetic parameter values for waste listed in Table 2 are selected based on the previous experiments on waste biodegradation     Fig. 3 presents the FS values for different friction angles calculated in the above two studies and by the proposed model without activating the gas flow and biodegradation submodels. The small differences among the three curves are due to the differences in the calculation method (LEM in Xu et al. (2012) and SRM in Feng et al. (2018) and this paper) and the modeling approach (FEM in Xu et al. (2012), FDM in Feng et al. (2018) and FVM in this paper). However, the differences are acceptable for engineering practice and the FS values calculated in this paper are generally conservative. Subsequently, the proposed model can then be applied in the following.

Leachate recirculation scenarios
The conceptual landfill in Fig. 2b is simulated in this part to investigate the effects of waste biodegradation and gas flow on the slope stability under leachate recirculation. Four scenarios are presented in Table 5, which can be achieved by activating the corresponding sub-model of the proposed model. To describe waste non-homogeneity with depth, the porosity, density and vertical intrinsic permeability of waste were assumed to linearly decrease from 0.5, 1500 kg/m 3 and 10 −12 m 2 at z = 50 m (top) to 0.3, 1800 kg/m 3 and 10 −14 m 2 at z = 0 m (bottom), respectively. A typical anisotropic coefficient of 5 was selected to consider waste anisotropy. Simulations with 11,400 finite-volume meshes were conducted on a workstation with 16 cores.

Gas pressure in landfill
Initially, the effect of waste biodegradation on the distribution of gas pressure under recirculation using horizontal trenches is investigated (P li = 49 kPa). Figure 4 compares the distributions of gas pressures at x = 60 m on the 5th, 10th, 20th, and 50th day in Scenarios 3 and 4. Overall, the gas pressures in both cases gradually increase over time due to the increased leachate saturation around the horizontal trench, and the gas pressures at the bottom of landfill are relatively higher since the leachate collection system is impermeable to gas. However, since the accelerated anaerobic degradation of waste will generate Dilation angle for waste ψ degree 0 Friction angle for waste φ degree 35 Xu et al. (2012) Waste cohesion c kPa 15 Fig. 3 Comparison of the relationships between FS and friction angle obtained by Feng et al. (2018), Xu et al. (2012) and the proposed model The proposed model in this study more gas, the gas pressure of Scenario 4 (0 ~ 16 kPa) is much larger than that of Scenario 3 (0 ~ 500 Pa) where waste biodegradation is neglected. Moreover, the gas pressure profiles of these two scenarios are totally different. In Scenario 3 without biodegradation, an evident turning point with a relatively high gas pressure can be observed at the front of infiltrated leachate. However, in Scenario 4, the turning point locates at the depth where the injected trench is installed because the higher leachate saturation around it can accelerate the gas generation and reduce the channel for gas migration at the same time (see Fig. 5). Figure 5 depicts the leachate pressure profiles (relative to atmospheric pressure) at x = 60 m on the 5th, 10th, 20th, and 50th day in Scenarios 1-4 and the saturation contour of scenario 4 is also included for comparison. The leachate pressure profiles of Scenario 1 (no gas flow and biodegradation) and Scenario 2 (no gas flow) are almost the same, indicating that the effect of generated leachate can be neglected when ignoring the gas flow. The difference between Scenario 1 and Scenario 3 (no biodegradation) is also small because no gas is generated when ignoring biodegradation. However, a significantly different leachate pressure profile can be observed in Scenario 4 (bio-hydro coupled) and the difference becomes more and more obvious over time. This demonstrates that the combination of waste biodegradation and gas flow obviously affects the redistribution of leachate during recirculation. In all the four scenarios, negative leachate pressures can be found at the areas near the top and bottom of the landfill, referring to an unsaturated zone due to leachate drainage. When the landfilled waste is saturated by the injected leachate, the leachate pressure becomes positive. On the 20 th day, the injected leachate from the buried horizontal trench can reach the bottom of the landfill in Scenarios 1-3, while in Scenario 4 there is still an unsaturated zone near the bottom since the generated gas impedes the leachate migration.

Leachate pressure in landfill
Many previous studies have neglected the effect of gas flow on leachate flow in an anaerobic bioreactor landfill. The results of this work manifest that this assumption is only reasonable in some cases where the organic content of landfilled waste is low and hence the waste biodegradation can be ignored. As a whole, both waste biodegradation and gas flow should be considered when analyzing the landfill gas and leachate flows, otherwise, the influence area and the performance of leachate recirculation would be overestimated.

Factor of safety of landfill slope
Before recirculation, the FS value of the modeled landfill slope is about 2.01, but with an injection pressure of 49 kPa, the FS value shows a reduction over recirculation time (Fig. 6), which is consistent with the increased pore pressures mentioned above. A small difference in the FS values of Scenarios 1-3 can be observed, corresponding to similar pore pressures presented in Fig. 5. Scenario 4 (bio-hydro coupled) has the smallest value of FS, and the FS difference between Scenarios 1 (no gas flow and biodegradation) and 4 gradually increases over time until reaching 0.35 on the 100th day (shadow region in Fig. 6). This difference is attributed to an increasing difference in the pore pressures (see Figs. 4 and 5) when considering the waste biodegradation and gas flow during leachate recirculation.
To further investigate the bio-hydro coupled effect on the landfill slope stability, Fig. 7 presents the FS value on the 100th day for 9 pairs of leachate injection pressure (P li = 49, 147, 245, and 343 kPa) and maximum anaerobic biodegradation rate (k n,max = 0, 0.01, 0.02, 0.03, 0.04, and 0.05 day −1 ). For the scenario with k n,max = 0 day −1 , the waste biodegradation is neglected. All the FS values calculated by the proposed model range between 1.31 and 1.96 in Fig. 7, within the normal range reported in the previous similar studies (Byun et al. 2019;Feng et al. 2018Feng et al. , 2020Xu et al. 2012). Both increasing the leachate injection pressure and maximum anaerobic biodegradation rate can reduce the FS value of a landfill slope maximally by 0.32 (at P li = 343 kPa) and 0.62 (at k n,max = 0.05 day −1 ), respectively, due to the increased pore leachate and gas pressures. Moreover, the FS difference between the bio-hydro coupled model and the   Fig. 7 Comparison between the FS values of different scenarios on the 100th day for 9 pairs of leachate injection pressure P li , and maximum anaerobic biodegradation rate k n,max model ignoring gas flow and biodegradation ranges between 0.33 and 0.37 (left sub-figure in Fig. 7). If only ignoring the gas flow, the waste biodegradation and the produced leachate have a very slight effect on the landfill slope stability, i.e., from 1.96 to 1.93 in the right subfigure of Fig. 7. However, the ignorance of both waste biodegradation and gas flow will largely overestimate the stability of a landfill slope by about 20-50% (e.g., (1.96-1.32)/1.32 = 46% for k N, max = 0.05 day −1 and P li = 49 kPa), especially when the landfilled waste has a high degradation rate and can then significantly increase the pore pressures.

Aeration scenarios
As one of the main strategies to accelerate waste biodegradation, aeration will increase the pore pressures of waste and then pose a threat to the landfill slope stability. In this part, the conceptual landfill in Fig. 2c is simulated using the proposed bio-hydro coupled model.
The porosity, density and vertical intrinsic permeability of waste were assumed to linearly decrease from 0.5, 1500 kg/m 3 and 10 −12 m 2 at z = 25 m (top) to 0.4, 1700 kg/ m 3 and 10 −13 m 2 at z = 0 m (bottom), respectively. A waste anisotropy coefficient of 5 was adopted. The aeration pressure P gi between 1 and 15 kPa was assumed base on engineering practice (Fytanidis and Voudrias 2014;Ritzkowski and Stegmann 2012;Townsend et al. 2015).
A leachate saturation of 0.5, an atmospheric gas pressure, and an oxygen fraction of 0% were assumed as the initial condition. Simulations with 80,000 finite-volume meshes were conducted on a workstation with 16 cores. Figure 8 shows the spatial distribution of gas pressure within the landfill in 3D view and cross-sectional view at y = 10 m on the 20th day of aeration with an injection pressure of 12 kPa. It can be observed that the gas pressure near the injection screen (about 9 ~ 12 kPa in red) and the bottom of the landfill (about 7 kPa in green) is significantly higher than that in other regions. The local accumulation of gas pressure within 3 m around the Fig. 8 The spatial distribution of gas pressure on the 20th day of aeration aeration well may be detrimental to the local stability of the landfilled waste. It should be noted that the gas pressure at the top was set as zero in this paper (in blue in Fig. 8) to represent an efficient gas collection system. Figure 9 shows the effect of aeration on the stability of a landfill slope in terms of the displacement at an unstable state. Compared to the scenario without air injection, an injection pressure of 12 kPa will slightly reduce the value of FS by 0.1 (4.5%) when the horizontal distance between the top of the injection screen and the slope surface (d) is 10 m. However, these two scenarios have different slope failure modes. For P gi = 0, the overall slope is unstable and the possible sliding surface develops from the left platform to the toe of the landfill slope. When P gi increases to 12 kPa, a local slope failure will occur before an overall instability mainly due to the accumulated gas pressures around the injection screen.
To further investigate the effect of aeration on the stability of a landfill slope, the changes of FS with air injection pressure P gi for different values of d (i.e., the horizontal distance between the top of the injection screen and the slope surface) are obtained and plotted in Fig. 10. As might be expected, the value of FS gradually decreases with increasing aeration pressure from 2.3 at P gi = 0 kPa (without aeration). The reduction in FS is slight at low aeration pressure, but becomes significant and then poses a threat to the slope stability when the aeration pressure exceeds a critical value, defined as the safe aeration pressure P gs below. The maximum reduction can reach up to 0.8 (about 53%, FS = 1.5) for P gi = 15 kPa and d = 5 m. Taking d = 7 m as an example, the FS remains stable around 2.3 and begins to show a significant decrease when P gi exceeds 3.8 kPa (red dash arrow in Fig. 10). Thus, P gs = 3.8 kPa when d = 7 m, that Possible slope failure mode in terms of the displacement at an unstable state: a P gi = 0 kPa; b P gi = 12 kPa is, an aeration pressure of P gi ≤ 3.8 kPa is recommended for d = 7 m during the operation of aeration in practice to benefit the landfill slope stability. By identifying the safe aeration pressure P gs for d = 5 ~ 10 m, a unique relationship between P gs and d can be observed in Fig. 10 (black solid line). P gs increases from 2 kPa at d = 5 m to 9.6 kPa at d = 10 m. These recommended aeration pressures regarding to the safety of the landfill slope are consistent with the values that are widely adopted in many aerobic projects, i.e., 2 ~ 8 kPa (Öncü et al. 2012;Raga and Cossu 2014;Ritzkowski et al. 2016;Townsend et al. 2015).
In summary, to maintain the stability of landfill slope under aeration, a smaller aeration pressure should be adopted for a shorter horizontal distance between the top of the injection screen and the slope surface, d. Based on the above numerical modeling, an aeration pressure exceeding 10 kPa is only allowed for d > 10 m. The design chart in Fig. 10 (shadow region) can provide some reference values for the air injection pressures at different values of d.

Conclusions
This study performs a three-dimensional numerical analysis for landfill slope stability during leachate recirculation and aeration using strength reduction method. The bio-hydro coupled processes of waste are simulated by a previously reported landfill coupled model programmed on the open-source platform OpenFOAM and then incorporated into the analysis. Compared with the previous studies, this study investigates the effects of waste biodegradation and leachate-gas flow on the stability of a landfill slope and provides some guidelines for leachate recirculation and aeration. The main findings are as follows: 1. The values of factor of safety calculated in this study using strength reduction method and finite volume method agree well with the previous related results, and the acceptable difference is due to the difference in the calculation method and modeling approach. 2. The landfill gas produced by waste biodegradation will significantly affect the values and the distribution of gas pressure in a landfill, and then affect the redistribution of leachate pressure and the slope stability during leachate recirculation. The ignorance of landfill gas flow during leachate recirculation is only reasonable when the organic content of landfilled waste is low and hence waste biodegradation can be ignored. Otherwise, the influence range and the performance of leachate recirculation would be overestimated. 3. Both increasing the leachate injection pressure and the maximum anaerobic biodegradation rate can reduce the FS value of the landfill slope maximally by 0.32 and 0.62, respectively, due to the increased pore leachate and gas pressures. If only ignoring the gas flow, the waste biodegradation and the produced leachate have a very slight effect on the landfill slope stability. However, the ignorance of both waste biodegradation and gas flow will significantly overestimate the stability of a landfill slope under leachate recirculation by about 20-50%, especially when the landfilled waste is easily degradable. 4. The FS of a landfill slope under aeration shows a significant reduction when the aeration pressure exceeds a critical value and this value is defined as safe aeration pressure. The maximum reduction can reach up to 0.8 (by about 53%). A relationship between the safe aeration pressure and the horizontal distance between the top of the injection screen and the slope surface, d, is proposed for preliminary design to avoid landfill slope failure during aeration. A smaller aeration pressure should be adopted when the distance d is small, and an aeration pressure exceeding 10 kPa is only allowable for d > 10 m.