Estimating Landfill Landslide Probability Using SAR Satellite Products: A Novel Approach

Abstract

This article presents a methodology for evaluating the susceptibility of landfill areas to develop landslides by analyzing Synthetic Aperture Radar (SAR) satellite products. The deformation velocity of the landfills is computed through the Persistent Scatterer Method on SAR imagery. These data, combined with a deformation model based on the shallow water equations (SWE), form the foundation for a Monte Carlo experiment that extrapolates the current state of the landfill into the future. The results of this simulation are then employed to determine the probability of a landslide occurrence. In order to validate the methodology effectiveness, a case study is conducted on a landfill in Zaldibar, Spain, revealing its effectiveness in estimating the probability of landfill landslides. This innovative approach emerges as an asset in large landfill management, acting as a proactive tool for identifying high-risk sites and preventing potential landslides, ultimately safeguarding human life and the environment. By providing insights into landslide probabilities, this study enhances decision-making processes and facilitates the development of intervention strategies in the domain of landfill risk assessment and management.

1. Introduction

The management of landfills poses a significant challenge in Europe, where an estimated 500,000 such sites exist, with a staggering 90% categorized as non-sanitary. Compounding the issue, approximately 80% of Europe’s landfills are publicly owned, making them susceptible to inadequate management practices [1]. These landfills, lacking essential environmental protection technologies, are on a trajectory towards demanding costly remediation efforts. Among the myriad concerns associated with these sites, landslides emerge as a paramount threat. Beyond the evident environmental toll, the uncontrolled conditions prevailing in these landfills have led to tragic incidents, exemplified by the 2020 Zaldibar disaster in Spain, resulting in two fatalities. The continent has witnessed more devastating examples, including the Lviv landslide in 2016 [2], and the 2023 Zagreb landslide. Outside Europe, incidents such as the Shenzhen landslide (2015) [3] or the Dona Juana landslide (1997) [4] serve as poignant examples, highlighting that this is a global problem with far-reaching consequences.

The profound environmental and human costs incurred by these events have spurred extensive research into their causes. Key factors identified in previous studies [5] include improper waste compaction practices, intense rainfall, and the accumulation of rainwater, which collectively heighten the susceptibility of municipal solid waste landfills to failure, often attributable to pore water pressure or insufficient compaction. This article delves into a novel approach for assessing and mitigating the risk of landfill landslides, offering a potential breakthrough in landfill management strategies.

Various models have been developed to study the dynamic behavior of landfills, aiming to predict potential landslides. The continuous deformation of landfills during and after the landfilling process, attributed to the high compressibility of waste, underscores the need for effective monitoring to mitigate the global repercussions of these failures. In a bid to estimate the safety factors of landfill slopes, Wang, Zhang, and Lin [6] proposed a two-dimensional limit analysis using the Morgenstern–Price method and EMU (energy method of upper bound theory). Their approach revolved around assessing the vertical profile composition of waste and its constituents. Building upon this, Zhang et al. [7] expanded the methodology through the incorporation of Monte Carlo simulations. This extension considered diverse factors such as varying water accumulation levels and geotechnical uncertainties within landfills. Zhang et al. calculated the probability distribution of the factor of safety (FOS).

In contrast, Xiu et al. [5] introduced a new approach to landslide study by treating landfill waste as a non-Newtonian fluid and employing Computational Fluid Dynamics (CFD) methods to analyze its evolution. This advanced methodology was successfully applied to investigate the 2015 Shenzhen “1220” landslide. Although the computational cost associated with this method surpasses that of conventional approaches, its notable advantage lies in the ability to achieve a significantly more accurate modeling of waste dynamics.

The accurate determination of landfill waste profiles is crucial when applying any of the mentioned methods. Traditionally, this involves on-site measurements or the utilization of UAV Lidar reconstructions [8]. However, this study suggests a departure from these conventional approaches by advocating for the use of cost-effective Synthetic Aperture Radar (SAR) satellite data as it was proposed by Papale et al. [9]. Unlike on-site measurements or UAV Lidar reconstructions, SAR data allow for the simultaneous monitoring of numerous landfills, presenting a more economically feasible and scalable solution for comprehensive risk assessment and management, at least in large landfills. The cost-effectiveness and efficiency of SAR satellite data make them well-suited for addressing the extensive nature of landfill sites across regions or countries, providing a practical alternative for researchers and practitioners in the field [10].

Given the limited resolution and precision of SAR data for adequately estimating the probability of landfill occurrences, a more robust approach involves simulating the dynamics of the landfill and conducting a statistical study through Monte Carlo experiments. To enhance the analysis with 3D effects while minimizing computational costs, the shallow water equations are employed [11,12]. This methodology has been previously utilized and validated in studies related to landslides. Specifically, the determination of the failure probability of the landfill takes into account uncertainties associated with shear strength parameters of the debris, encompassing considerations of the friction angle and cohesion.

The paper’s structure is outlined as follows. Section 2 provides an explanation of the materials and methods employed in this study. This includes a description of the satellite data used (Section 2.1) and a comprehensive overview of the algorithm (Section 2.2), with detailed explanations of the main relevant methods provided in each subsection. The verification and validation of the landfill deformation propagator are conducted in Section 3.1 and Section 3.2, respectively, while the algorithm is tested against the real case of the Zaldibar landslide in Section 3.3. Lastly, Section 4 presents a detailed discussion of the main conclusions drawn from the study.

2. Materials and Methods

2.1. Satellite Data

When choosing the most suitable SAR satellite for monitoring purposes, critical factors to consider include spatial resolution, flexible data access, and revisit time.

Table 1. Selected SAR satellites.

	Frequency Range	Resolution	Frequency
Sentinel 1	C-Band (5.4 GHz)	5 × 20 m	12 days 1
COSMO SkyMed	X-Band (9.65 GHz)	1 m	5 days
RADARSAT	C-Band (5.4 GHz)	1 × 3 m	24 days
Terra SAR X	X-Band (9.65 GHz)	1–16 m	11 days
ALOS 2	L-Band (1.24 GHz)	3 m	14 days
RISAT	X-Band (9.65 GHz)	1 m	4 days

¹ Before the Sentinel-1B breakdown (December 2021), its revisit frequency was every 6 days.

presents a concise overview of key SAR satellites, outlining their frequency range, resolution, and data acquisition frequency. In the context of this scenario, all available satellites demonstrate resolutions ranging from 1 to 20 m, a range proven sufficient for landfill monitoring in numerous previous studies [7]. Emphasizing the other key parameters, namely revisit time and data availability, Sentinel-1 emerges as the optimal choice for the current study. Notably, Sentinel-1 has consistently proven to be the preferred solution in the most recent studies focused on monitoring surface deformation in landfills [10].

For this study, 13 Sentinel-1 Single Look Complex (SLC) images acquired in Interferometric Wide Swath (IW) mode, sub-swath 1, with VV polarization have been employed. The dataset spans the period from 1 December 2018 to 13 January 2020.

2.2. Proposed Methodology

The algorithm for predicting landslides in landfills through SAR image analysis is depicted in Figure 1, outlining a methodology that comprises four key phases. To initiate the process, SAR satellite data for the landfill site are acquired. The preferred method for accessing Sentinel-1 products is through the Copernicus Open Access Hub, a web-based portal designed for data retrieval from the Copernicus Sentinel missions. The obtained data are then downloaded in SAFE format.

Subsequently, the algorithm progresses to data preprocessing, a stage dedicated to refining the raw SAR data to make them suitable for in-depth analysis. Tasks such as radiometric calibration, speckle filtering, and image registration are carried out during this phase. The chosen technique, the multitemporal InSAR method known as Persistent Scatterer Interferometry (PSI), is employed to address the limitations of conventional InSAR methods [13], thereby enhancing accuracy in detecting ground deformation.

With the ultimate goal of determining landslide probability, a Monte Carlo experiment is conducted. This involves varying the dynamic properties of various debris materials [14] and projecting the current state of the landfill into the future. Previous studies, such as Dai et al. [15], have identified the unit weight and moisture content as the most significant properties for the dynamics of debris deformation. The unit weight varies from 6 kN/m³ near the surface to approximately 13 kN/m³ at depths exceeding 45 m. The typical moisture content ranges between 0.75 and 0.95. Some combinations of these properties may lead to a landslide if the landfill is deemed unstable, based on the initial state of the landfill.

Due to the limitations of InSAR techniques in detecting significant changes in landfill volumes, such as when new waste deposits are made, determining the current state of the landfill—specifically, computing the exact amount of Municipal Solid Waste (MSW) contained—poses a challenge. Various options exist, including on-site measurements, historical data analysis, UAV surveys, and more. However, in this study, the chosen approach involves estimating the initial conditions through the resolution of the inverse problem (refer to Figure 2). By utilizing the known current deformation velocities, the volume of waste responsible for generating those velocities is estimated. This estimation is achieved by applying a nonlinear minimization solver to the equations of the shallow water equations (SWE) model.

Upon completion of the Monte Carlo experiment and the availability of results, a postprocessing step is carried out to identify the most critical areas susceptible to landslides. This involves analyzing the number of realizations in which the debris flow reaches a particular region, aiding in the detection of areas near the landslide that are at risk. It is essential to note that the quality of the results hinges on several parameters, including the spatial resolution of SAR measurements, the proper design of uncertainty variables employed in the Monte Carlo experiment, and the utilization of an accurate digital terrain model for the area of interest. A meticulous consideration of these factors contributes to the reliability and precision of the findings, enhancing the effectiveness of the analysis in identifying potential landslide-prone areas.

2.2.1. Persistent Scatterer Interferometry (PSI) Analysis

The prescribed working procedure, based on the PSI technique, is described by Nam et al. [16]. In PSI, the first step involves selecting a master image from a set of available (n + 1) images. Interferograms are then generated by comparing the master image with the $n$ slave images, which have been acquired on different dates. Each interferogram is characterized by a specific perpendicular baseline, representing the distance between the satellite positions. By accurately determining the satellite’s position on the acquisition dates, it becomes possible to remove components related to the perpendicular baseline, such as flat earth and topographic effects, from the interferogram phase.

The initial step involves processing the data using ‘snap2stamps’ (v1.01) [17,18] a software that facilitates pre-processing through the SNAP GPT (Graphic Processing Tool, v9.0.0). This processing prepares the data by converting them to the gamma format, which is compatible with Stamps for the application of the PSI method. The following steps are applied to the Sentinel-1 products, following the StaMPS workflow. More information about each step can be found in the original work [19]:

• Master/slave selection and splitting.
• Co-registration and Interferogram computation: This step conducts co-registration between the master image and each of the slave images in sequence. The interferograms corresponding to each pair of master and slave images are then generated prior to removing the flat-earth phase, i.e., the phase associated with the ellipsoid. Debursting for both the SLC images and the differential interferograms is applied to remove horizontal stripes. The topographical phase is then simulated using the 3 s Shuttle Radar Topography Mission (SRTM) Digital Terrain Model (DTM), which is downloaded automatically by SNAP. This topographical component is then removed from the interferograms. In this step, a subset of the previously selected research area is chosen.
• StaMPS export for PSI. The export step using StaMPS is performed using the operator with the same name, and the inputs required for this step are: (i) the co-registered master–slave pair, (ii) the corresponding interferogram with the elevation and latitude and longitude bands that have been orthorectified.
• Data load.
• Estimate phase noise.
• PS selection: Pixels are selected on the basis of their noise characteristics; this step also estimates the percentage of random (non-PS) pixels in a scene from which the density per km² can be obtained.
• PS weeding: In the previous step, pixels are selected and filtered, discarding those that are caused by contributions from neighboring ground resolution elements and those considered too noisy. Data for the selected pixels are stored in new workspaces.
• Phase correction: The wrapped phase of the selected pixels is corrected for spatially uncorrelated look angle (DEM) error. At the end of this step, the patches are merged.
• Phase unwrapping.
• Estimate spatially uncorrelated look angle error: spatially uncorrelated look angle (SULA) error was calculated in Step 3 and removed in Step 5. In Step 7, spatially correlated look angle (SCLA) error is calculated which is due almost exclusively to spatially correlated DEM error (this includes error in the DEM itself, and incorrect mapping of the DEM into radar coordinates). Master atmosphere and orbit error (AOE) phase are estimated simultaneously.
• Atmospheric filtering.

Following these steps, the deformation extraction process analyzes the unwrapped interferograms to calculate the surface deformation that has occurred between SAR image acquisitions.

SAR velocities are inherently measured along the Line of Sight (LOS) direction. To gain a more comprehensive understanding of surface movement, it is necessary to transform them into local geodetic systems (e.g., north, east, up). This study employs the method described in [20] for this transformation.

However, due to Sentinel-1’s acquisition geometry, calculating horizontal velocity components in all directions (e.g., north–south) is limited. Landfill landslides are primarily driven by factors like waste compression and water content variations. Therefore, this study initially focuses on the vertical displacement component, which is readily derived from the angle of incidence and satellite orbit. This approach prioritizes the dominant movement component during the initial landslide phase, neglecting the potential contribution of horizontal components.

2.2.2. Deformation Model

To assess landslide likelihood in potentially hazardous areas (identified using the InSAR technique), a 2D shallow water equation solver is implemented in Matlab, utilizing the finite volume method and a total-variation-diminishing (TVD) scheme [21]. This model serves two purposes: to predict landslide probability by simulating future landfill states and to estimate the landfill’s initial waste volume through an iterative method incorporating the deformation model and deformation velocities calculated in Section 2.2 and Section 3.3.

The characteristics of mass movements in landslides can be approximated by the Reynolds-averaged Navier–Stokes equations, expressed in the differential conservation form as:

$$\mathrm{Mass} \mathrm{conservation}: \frac{\partial \rho }{\partial t}+\nabla \cdot \rho \tilde{\mathrm{u}}=0$$

(1)

$$\mathrm{Momentum} \mathrm{conservation}: \frac{\partial \rho }{\partial t}+\nabla \cdot \rho \mathit{uu}=\nabla \cdot \mathsf{\tau }+\rho g$$

(2)

where $\rho$ is mass density, $g$ is gravitational acceleration, $t$ is the time, $u$ is the mass velocity vector (the variables are averaged in the vertical direction) and $\mathsf{\tau }$ is the stress tensor.

To solve Equations (1) and (2), they can be reduced to a two-dimensional shallow water form by integrating over the vertical axis. This simplification is valid if the movement of the mass in the horizontal direction outweighs that of the vertical direction [22,23]. It results in a hyperbolic partial differential equation whose conservative differential form is:

$$\frac{\partial \mathrm{Q}}{\partial t}+\frac{\partial \mathrm{F}}{\partial x}+\frac{\partial \mathrm{G}}{\partial y}=\mathrm{S}$$

(3)

in which:

$$\begin{array}{c}\mathrm{Q}=\left[\begin{array}{c}h\\ hu\\ hv\end{array}\right] \mathrm{F}=\left[\begin{array}{c}hu\\ h{u}^2+\frac{g{h}^2}{2}\\ huv\end{array}\right]\\ \mathrm{G}=\left[\begin{array}{c}hv\\ huv\\ h{v}^2+\frac{g{h}^2}{2}\end{array}\right] \mathrm{S}=\left[\begin{array}{c}0\\ s_{0x}+s_{fx}\\ s_{0y}+s_{fy}\end{array}\right]\end{array}$$

where $\mathrm{Q}$ is the state vector and holds both the terrain depth $h$, and the depth-averaged velocity cartesian components $u$ and $v$; $\mathrm{F}$ and $\mathrm{G}$ are the flux vectors in the cartesian components; and $\mathrm{S}$ is a source term. For the latter matrix, the formulation found in [24] is used.

There are two parameters exclusive to the source term. Firstly, the bed slopes, $s_{0x}$ and $s_{0y}$, that depend on the local slopes ($\frac{\partial z_f}{\partial x}$, $\frac{\partial z_f}{\partial y}$) of the ground:

$$s_{0x}=-\frac{gh\frac{\partial z_f}{\partial x}}{{\left(1+{\frac{\partial z_f}{\partial x}}^2+{\frac{\partial z_f}{\partial y}}^2\right)}^{\frac{1}{2}}} s_{0y}=-\frac{gh \frac{\partial z_f}{\partial y}}{{\left(1+{\frac{\partial z_f}{\partial x}}^2+{\frac{\partial z_f}{\partial y}}^2\right)}^{\frac{1}{2}}}$$

Secondly, the bed friction slopes, $s_{fx}$ and $s_{fy}$, which are expressed as a function of depth, velocity and, most important, the effective friction coefficient $\mu$ following the Coulomb friction model [25]:

$$s_{fx}=-\frac{\mu gh}{{\left(1+{\frac{\partial z_f}{\partial x}}^2+{\frac{\partial z_f}{\partial y}}^2\right)}^{\frac{1}{2}}}\cdot \frac{u}{\sqrt{u^2+v^2}} s_{fy}=-\frac{\mu gh}{{\left(1+{\frac{\partial z_f}{\partial x}}^2+{\frac{\partial z_f}{\partial y}}^2\right)}^{\frac{1}{2}}}\cdot \frac{v}{\sqrt{u^2+v^2}}$$

Numerical Implementation

The implementation for the numerical simulation of the SWE follows the procedure described in [21]. This model has been chosen because it manages to maintain conservative properties as well as solving discontinuities without spurious oscillations while not introducing too much numerical dissipation.

Applying the divergence theorem over a control volume $\mathsf{\Omega }$ to Equation (1), it can be expressed as

$${{\displaystyle \iint }}_{\mathsf{\Omega }}\frac{\partial \mathrm{Q}}{\partial t}\mathrm{d}W+{{\displaystyle \int }}_{\partial \mathsf{\Omega }}\mathrm{E}·\mathrm{n} \mathrm{d}l={{\displaystyle \iint }}_{\mathsf{\Omega }}\mathrm{S} \mathrm{d}W,$$

(4)

in which $\mathrm{n}$ is the outward unit vector normal to the boundary $\partial \mathsf{\Omega }$; $\mathrm{d}W$ and $\mathrm{d}l$ are the area and arc elements; and the integrand $\mathrm{E}·\mathrm{n}$ is the normal flux vector in which $\mathrm{E}={\left[\mathrm{F}, \mathrm{G}\right]}^T$.

Assuming that $\mathrm{Q}$ is constant all over the cell and represents the average value at the center, the discretization in the Finite Volume Method is:

$$A\frac{dQ}{dt}+{\displaystyle \sum }_{m=1}^M{\mathrm{E}}_{\mathrm{n}}^m{L}^m=A·\mathrm{S}$$

(5)

where $A$ is the area of the cell, $m$ is the index for the side, $M$ is the total number of sides for a cell, ${\mathrm{E}}_{\mathrm{n}}^m$ is the intercell flux normal to the side, and $L^m$ is the length of the side. Applying the rotational invariance [19], the intercell flux normal to each side is defined as

$${\mathrm{E}}_{\mathrm{n}}\left(\mathrm{Q}\right)=\mathrm{F}\cos \mathsf{\Phi }+\mathrm{G}\sin \varphi =\mathrm{T}{\left(\mathsf{\Phi }\right)}^{-1}\mathrm{F}\left[\overline{\mathrm{Q}}\right]=\mathrm{T}{\left(\mathsf{\Phi }\right)}^{-1}\mathrm{F}\left(\overline{\mathrm{Q}}\right)$$

(6)

where $\mathsf{\Phi }$ is the angle measured counterclockwise from the $x$ axis from the vector $\mathrm{n}$, $\overline{\mathrm{Q}}=\mathrm{T}\left(\mathsf{\Phi }\right)\mathrm{Q}$ is the vector with the transformed variables to the coordinates $\overline{\mathrm{Q}}={\left(h, h{u}_n,h{v}_t\right)}^{\mathrm{T}}$, being:

$$\left\{\begin{array}{c}{u}_n=u\cos \mathsf{\Phi }+v\sin \mathsf{\Phi }\\ v_t=v\cos \mathsf{\Phi }-u\sin \mathsf{\Phi }\end{array}\right.$$

(7)

in which $u_n$ is perpendicular to the side and $v_t$ is parallel. Substituting Equation (6) into Equation (5):

$$A\frac{\mathrm{dQ}}{\mathrm{d}t}+{\displaystyle \sum }_{m=1}^M\mathrm{T}{\left(\mathsf{\Phi }\right)}^{-1}\mathrm{F}\left(\overline{\mathrm{Q}}\right)L^m=A·\mathrm{S}$$

(8)

Once discretized, the objective of the TDV scheme is to estimate the flux between every cell interface, $\mathrm{F}\left(\overline{\mathrm{Q}}\right)$, with two similar steps to preserve second-order accuracy in both time and space.

$${\hat{\mathrm{Q}}}_{i,j}={\mathrm{Q}}_{i,j}^n-\frac{\mathsf{\Delta }t}{A}{\left[{\displaystyle \sum }_{m=1}^M\mathrm{T}{\left(\mathsf{\Phi }\right)}^{-1}·{\mathrm{F}}^{\left(1\right)}\left(\overline{Q}\right)·L^m \right]}_{i,j}+\mathsf{\Delta }t·{\mathrm{S}}_{i,j}^n.$$

(9)

The first step predicts the intercell flow by taking finite differences of Equation (6). Note that the index $i,j$ is used to refer each cell individually and $n$ to denote the time index. The new term ${\mathrm{F}}^{\left(1\right)}\left(\overline{\mathrm{Q}}\right)$ refers to the first-order numerical flux. It is different in each side of the cell. In $(i+1/2,j$), which corresponds to the interface between $\left(i,j\right)$ and $\left(i+1,j\right)$, it can be estimated as

$${\mathrm{F}}^{\left(1\right)}\left(\overline{\mathrm{Q}}\right)={\mathrm{F}}_{\mathrm{LR}}={\mathrm{F}}_{\mathrm{L}}^{+}+{\mathrm{F}}_{\mathrm{R}}^{-}={\mathrm{F}}_{i,j}^{+}+{\mathrm{F}}_{i+1,j}^{-}.$$

(10)

Here, we are decomposing the outward flow into positive and negative parts with a technique known as van Leer splitting (denoted as VLS).

$${\mathrm{Q}}_{i,j}^{n+1}={\mathrm{Q}}_{i,j}^n-\frac{\mathsf{\Delta }t}{A}{\left[{\displaystyle \sum }_{m=1}^M\mathrm{T}{\left(\mathsf{\Phi }\right)}^{-1}·{\mathrm{F}}^{\left(2\right)}\left(\overline{Q}\right)·L^m \right]}_{i,j}+\mathsf{\Delta }t·{\hat{\mathrm{S}}}_{i,j}.$$

(11)

In the corrector step, the value for the matrix ${\mathrm{Q}}_{i,j}$ in the next time step $n$ is computed. The procedure is very similar to the predictor with two exceptions: (1) the source term matrix ${\mathrm{S}}_{i,j}$ is reformulated with the predicted values for $h$, $u$ and $v$; and (2) a second-order numerical flux, ${\mathrm{F}}^{\left(2\right)}\left(\overline{\mathrm{Q}}\right)$ is estimated. For this flux, an antidiffusive term has been added to the previous first-order numerical flux ${\mathrm{F}}_{\mathrm{LR}}$.

$${\mathrm{F}}^{\left(2\right)}\left(\overline{\mathrm{Q}}\right)={\mathrm{F}}_{\mathrm{LR}}+\frac{1}{2}\left[\varphi \left(r_{i+1/2,j}^{+}\right)w_{i+1/2,j}^{+}-\varphi \left(r_{i+1/2,j}^{-}\right)w_{i+1/2,j}^{-}\right],$$

(12)

with

$$r_{i+1/2,j}^{+}=\frac{w_{i-1/2,j}^{+}}{w_{i+1/2,j}^{+}}, r_{i+1/2,j}^{-}=\frac{w_{i+3/2,j}^{-}}{w_{i+1/2,j}^{-}},$$

$$w_{i+1/2,j}^{+}={\hat{\mathrm{F}}}_{i+1,j}^{+}-{\mathrm{F}}_{i,j}^{+}, w_{i+1/2,j}^{-}={\mathrm{F}}_{i+1,j}^{+}-{\hat{\mathrm{F}}}_{i,j}^{+}, w_{i+3/2,j}^{-}={\mathrm{F}}_{i+2,j}^{+}-{\hat{\mathrm{F}}}_{i+1,j}^{+},$$

$$\varphi \left(r\right)=\frac{r+\left|r\right|}{1+\left|r\right|},$$

where $\hat{\mathrm{F}}$ stands for the matrix $\mathrm{F}$ constructed with the predicted values, and $\varphi \left(r\right)$ is known as a van Leer limiter and is added to avoid unwanted oscillations propagating in further iterations.

2.2.3. Monte Carlo Experiment Design

The propagation of deformation velocities in landfills is subject to various uncertainties, primarily arising from factors such as the low resolution of SAR images, inaccuracies in the digital terrain model, errors in the propagation of deformation models, and the inherent uncertainty about the properties of the landfill’s waste material.

The first source of uncertainty is the low resolution of SAR images, which can limit the level of detail captured in the measurements of deformation velocities. This can result in reduced accuracy when detecting and quantifying subtle deformations within the landfill. Inaccuracies in the digital terrain model used for analyzing the landfills can also contribute to uncertainties in the propagation of deformation velocities. The digital terrain model serves as a representation of the land’s surface, and any errors or inaccuracies in this model can impact the analysis results and subsequent deformation predictions.

Additionally, errors in the propagation of deformation models introduce further uncertainties. The SWE model is utilized to interpret the collected data and predict the behavior of the landfill. However, this model is based on assumptions and simplifications, which can lead to propagation errors during the analysis process. Lastly, the inherent uncertainty about the properties of the landfill’s waste material can affect the propagation of deformation velocities. Municipal solid waste consists of diverse components with varying properties, such as density, moisture content, and mechanical behavior. Characterizing these properties accurately is challenging, and the uncertainty surrounding them can introduce uncertainties in the analysis and prediction of deformation velocities. Considering that incorporating all uncertainties is challenging, a proposed approach is to model the uncertainties associated with the material properties by incorporating a sufficient security margin that encompasses all other sources of uncertainty.

By modeling the uncertainties in the material properties, such as the mechanical behavior, density, and moisture content of the landfill, a safety margin can be included to account for other factors contributing to the overall uncertainty. This approach acknowledges that uncertainties stemming from the low resolution of SAR images, inaccuracies in the digital terrain model, and errors in the deformation model propagation are difficult to quantify precisely. Therefore, by incorporating a security margin based on the uncertainties about the material properties, a more comprehensive assessment can be achieved.

The study conducted by Dai et al. [15] emphasized the importance of considering the unit weight and water content of MSW in landfill deformation analyses. Kavazanjian et al. [26] provided widely referenced MSW unit weight profiles, indicating that the unit weight ranges from around 6 kN/m³ near the surface to approximately 13 kN/m³ at depths exceeding 45 m.

On the other hand, the work by Hossain [27] provided valuable insights into the water content variations within different types of MSW. High water content has been identified as a significant contributor to flow failures in landfills, as evidenced by the findings at the Payatas Landfill [28] and it is related to the shear strength distribution. To study the flow behavior of MSW, the research focused on preparing MSW simulants with elevated water content. Moisture contents ranging from 0.75 to 0.95 were selected to create MSW simulants with higher mobility. This enabled the observation and analysis of the flow behavior under these specific conditions.

As specified in the deformation model submodule, two key parameters that characterize the material are the equivalent viscosity coefficient and the density. For the Monte Carlo experiment, considering the comprehensive studies conducted earlier, the density will range from 600 kg/m³ to 1300 kg/m³. Additionally, the effect of water content will be accounted for in the equivalent viscosity coefficient, encompassing the range of 0.75–0.95.

3. Results

3.1. Verification of the Solver

The initial verification of the code has been performed by applying it to a standard test case involving the partial breach of a dam. The primary objective of this validation was to ensure that the implemented equations are correctly implemented within the code. Several researchers [29,30,31,32] have reported on this test case. The computational domain, depicted in Figure 3, consists of a 200 m × 200 m closed area with non-slip walls, featuring a 75 m wide breach and a 10 m thick dam structure in the direction of the flow. The initial conditions involve a stagnant water body with a height of 10 m in the left half and 5 m in the right half, including the section representing the dam/sluice gate.

The computation is performed for a duration of 7.2 s, at which point a well-developed bore is present in the central portion, and the wavefront has reached one bank of the channel. Figure 4 illustrates the three-dimensional water surface at the final time and displays iso-value lines of the water height in our calculations (ULE) compared to the aforementioned authors.

Additionally, Figure 5 showcases the computed water surface along cuts at x = 110 m and y = 130 m, with numerical results from [29,32,33,34] also provided for comparison. Fennema employed a shallow water model solved using an implicit finite difference method, while Biscarini [33] also implemented the full 3D Navier–Stokes multiphase equations in OpenFOAM. The overall shape of the water height demonstrates good agreement with the reference results.

3.2. Validation of the Capability to Predict the Landfill Deformation

The objective of the upcoming validation is to establish the accuracy of these equations in effectively depicting landfill deformation. To accomplish this, the study utilizes experimental data from [15]. In their research, the authors conducted a series of model tests to examine the dynamic behavior of municipal solid waste that experienced collapses resulting from landfill slope failures. These tests involved performing ring shear tests on a simulant of MSW under varying shear rates.

The inclined base angle of the model box was adjusted between 0 deg and 15 deg. Once the baffle is lifted, the simulated MSW experiences instantaneous collapse, flows down along the inclined base, and eventually settles along the bottom of the model box.

In order to carry out the validation, the initial experimental conditions were replicated to compare the distance traveled by the MRW for different slope angles. This comparison aimed to assess the accuracy and consistency of the model’s predictions with respect to the observed behavior in the experiments. By utilizing the test case with a slope of 0 deg, the equivalent viscosity of the model was calibrated accordingly. Following this adjustment, the slope of the model was progressively increased, and a comprehensive comparison was made against the experimental data. The outcomes demonstrated a commendable performance of the numerical SWE model, as shown in

Table 2. Runaway distance predicted by the 2D shallow water equations against the model test results [15,35] model.

Slope α (deg)	Model Test (cm)	Huang and Cheng (cm)	Error (%)	2D Shallow Water (cm)	Error (%)
0	53.5	58.27	8.92	53.5	0.00
5	60.8	64.14	5.49	62.6	2.81
10	76.5	69.93	−8.59	75.76	−1.06
15	80.2	75.41	−5.97	82.8	3.45

.

3.3. Case of Study: Zaldibar Landslide

Zaldibar is a municipality in the province of Vizcaya, Spain. On 6 February 2020, a massive landslide occurred at a landfill located on the slope of Mount Egoarbitza. The landfill was managed by Verter Recycling 2002 and was used for the disposal of industrial and construction waste.

The landslide buried two workers under tons of waste [36]. The cause of the landslide was attributed to the accumulation of waste and improper management practices. The waste had been inadequately compacted and had been accumulating for years, leading to unstable conditions. The incident triggered a major rescue operation to locate and recover the missing workers. Unfortunately, it took weeks to retrieve their bodies due to the hazardous conditions and the need for careful stabilization of the site. The Zaldibar landfill landslide also had severe environmental consequences. The collapse resulted in the release of a large amount of waste, including hazardous materials, into the nearby river, causing significant pollution and environmental damage.

The collapse caused the sliding of 800,000 cubic meters of debris down the entire side of the mountain. This, in turn, caused a large landslide 160 m wide and 330 m long, for a total area of 2780 square meters. This can be seen graphically in Figure 6. The waste spilled into the nearby river Ibaizabal, causing severe environmental contamination, and posing risks to the ecosystem and human health. The environmental impact of the incident led to the temporary closure of fishing and shellfish harvesting in the affected area, as well as the deployment of containment booms and other measures to mitigate the spread of pollutants. Cleanup efforts were initiated to remove the waste and restore the affected ecosystems.

In the aftermath of the incident, investigations were conducted to determine the causes of the landslide and the responsible parties. It was found that the accumulation of waste at the landfill had exceeded its designed capacity, leading to unstable conditions. Additionally, inadequate waste management practices, including improper waste compaction and lack of monitoring, were identified as contributing factors.

The Zaldibar landslide is a reminder of the risks associated with poorly managed landfills. It is important to ensure that landfills are properly designed, constructed, and operated to minimize the risk of landslides and other environmental hazards.

Following the application of the PSI technique, vertical displacement values for each of the Persistent Scatterers are determined, as depicted in Figure 7. Persistent Scatterers are typically identified in regions characterized by stable and coherent radar reflectors that exhibit minimal changes over time. These areas encompass various environments, such as urban settings with consistent building structures, rocky terrains with unchanging surfaces, agricultural fields featuring stable vegetation, desert landscapes with stationary sand or rock features, coastal regions hosting fixed coastal infrastructure, and industrial sites with steady installations [37].

Figure 6. Schematic of the Zaldibar landslide. (A) Landfill, (B) old quarry, (C) maximum extent of the landslide, (D) fires caused by methane. Adapted from [38].

Figure 6. Schematic of the Zaldibar landslide. (A) Landfill, (B) old quarry, (C) maximum extent of the landslide, (D) fires caused by methane. Adapted from [38].

While the focus of our paper centers on the outcomes related to the landfill for clarity, it is important to acknowledge the wider scope of our analysis. The processing area extended across the geographic coordinates between longitude (−2.75°, −2.25°) and latitude (43°, 42.3°). Within this expansive range, our detailed Persistent Scatterer Interferometry (PSI) analysis revealed a mean PS point density of 19 points per square kilometer.

Figure 8 depicts the isocontours of the interpolated deformation velocity around the Persistent Scatterers. In this figure, it is evident that only one point within the landfill area approaches a neutral state, while the remainder exhibits subsidence. The most substantial variations are observed in the southeast and northwest regions, corresponding to the data collected for mesh generation. This information holds significant value in enhancing our comprehension of the landfill area’s dynamics and can inform decision-making processes related to its management.

The presence of a single point approaching a neutral state suggests a potential collapse within the landfill area. This collapse could be attributed to various factors, including waste accumulation or the presence of materials with distinct properties. The observed variations in the southeast and northwest points indicate ground movements that have the potential to impact the stability of the landfill. Consequently, it is crucial to maintain ongoing monitoring of these areas to assess associated risks and facilitate well-informed decision-making processes.

Due to the limited density of Persistent Scatterers within landfill, drawing meaningful conclusions about their stability is challenging. As mentioned before, as an alternative approach, the focus was shifted to propagate the current volume of debris within the landfill into the future. This volume can, in fact, be computed from the deformation velocities previously measured.

Initially, the landfill area is delineated by employing a quadrilateral shape (see Figure 9), although more intricate geometries may be considered in the future. Subsequently, an initial estimate of the waste volume within these boundaries is made. This initial volume undergoes adjustment through the application of a third-order polynomial. These initial conditions are then propagated using the shallow water conditions and a digital terrain model. If the deformation velocities do not align closely with the measurements, the waste volume is adjusted iteratively until convergence is achieved. The nonlinear optimization problem is solved applying the standard Nelder–Mead simplex direct search method [39]. This methodology enabled the estimation of the landfill’s initial volume (the outcome for the Zaldibar landfill is depicted in Figure 10), offering crucial insights into its characteristics.

Following this, a Monte Carlo experiment comprising 50 realizations was developed. Within each realization, variations in waste density and water content were introduced, guided by data from Dai, Huang, and Jiang [15]. The simulations continued until the velocity deformations reached a negligible level.

Challenges with spatial resolution have been encountered that have affected the quality of the results. It appears that the algorithm could benefit significantly from higher resolution data, which would enhance our ability to model the initial volume of debris more accurately. Figure 11 displays the mean waste height deviation resulting from the Monte Carlo experiment, offering insights into variations across the landfill. Meanwhile, Figure 12, positioned on the right, illustrates the standard deviation. Notably, there are regions where the standard deviation markedly deviates from the norm observed in the rest of the landfill area. These disparities suggest the presence of localized irregularities or unique characteristics that warrant further investigation and analysis.

To assess the relevance of the results in predicting landslides, it is essential to compare them with the actual landslide event. The landslide occurred in the northern part of the landfill, with the debris flowing northward until reaching the highway. As depicted in Figure 13, the Monte Carlo method predicts that the most likely scenario is the debris flowing in the same direction as the actual landslide, with an approximate 50% probability of reaching the highway. What is particularly intriguing is that, among all the possible paths the debris could take while descending, our method accurately identifies the real trajectory of the debris, avoiding other downward directions. This highlights the effectiveness of our approach in replicating the observed landslide dynamics.

Nevertheless, it is probable that the method is currently overestimating the likelihood of a landslide occurrence. In all the experiment realizations, debris flow is consistently detected. This discrepancy may be rectified in the future by incorporating more realistic values for debris properties and enhancing the modeling of the initial waste volume. These refinements will contribute to a more accurate representation of landslide probabilities in subsequent analyses. To improve the accuracy of our predictive methods, it is essential to analyze additional landfill landslides.

4. Discussion and Conclusions

Improper management of landfills can result in landslides, causing significant environmental problems such as soil erosion, water contamination, and habitat destruction. Additionally, these landslides pose a serious threat to human life and may lead to fatalities. Given that a substantial proportion of these landslides occur in illegal or poorly managed landfills, there is a pressing need for an economical method by which governmental agencies can effectively control a large number of uncollaborative waste disposal sites. Addressing this issue is crucial for mitigating environmental risks and safeguarding public safety.

In this context, the proposed method boasts several advantages over preceding approaches. As it only requires SAR satellite data, it can be applied across diverse locations on Earth with a low cost, eliminating the need for collaboration with landfill managers. The deformation propagation method, validated and proven effective, imposes minimal computational burdens. Furthermore, it demonstrates the capability to integrate the 3D effects of the digital terrain model alongside key properties of the waste.

Sentinel-1 SAR products were selected for their applicability in this study, with SAR images of the Zaldibar landfill successfully acquired using the “Copernicus Open Access Hub”. The chosen methodology involves processing SAR images with the PSI algorithm, known for its precision and sensitivity to ground deformation.

In the application of this methodology to the Zaldibar landfill, a comprehensive validation was conducted to assess its effectiveness in landslide identification. A notable challenge arose from the low density of Persistent Scatterers (PS) within the landfill, particularly along its perimeter. This limitation, influenced by debris behavior in relation to radar measurements, posed constraints on the methods employed.

The algorithm detected a high probability of collapse and a substantial hazard of impacting the highway, which could have been used to prevent fatalities and environmental damage, such as pollution, toxic spills, and water contamination. The Monte Carlo method calculated the most probable scenario of the debris moving along the same trajectory as the actual landslide, with an approximate 50% probability of reaching the highway.

While future analyses should consider additional landfill landslide cases to enhance predictive accuracy and refine methods for initial waste volume estimation, as well as explore the impact of spatial resolution, this study offers valuable insights into the potential and limitations of SAR technology for monitoring landslides in complex environmental conditions. The findings contribute to the ongoing efforts to advance the understanding and application of SAR technology in the context of landfill monitoring, laying the groundwork for further improvements and optimizations in predictive modeling and analysis.