Solving geophysical flow problems in Luxembourg: SPH constitutive review

Abstract

In this research work, a constitutive mathematical review has been carried out to study the application of smoothed-particle hydrodynamics (SPH) in the modeling of geophysical flows like landslide, debris flows and stability failure problems across the world with particular focus on the landslide-associated geohazards happening in the southern region of Luxembourg. A particular interest is drawn to this region due to the research position taken by Stefan Van Baars suggesting the sufficiency of using friction angle of the geophysical flow-prone southern region to model landslide problems. This position invariably suggests that to model geophysical flows in any other region faced with this problem, friction angle data should have been enough. Meanwhile the present study has presented that climate factors and soil behavior variables beyond friction angle are required to solve geophysical flow problems and proposed the mathematical SPH technique as being the most efficient method to model geophysical flow problems due to its versatility and flexibility. Also, the (SPH) method can be used to simulate slope movements with large deformation and requires comprehensive profiling of the environmental and field data of the studied region. Generally, the present research has shown from previous efforts made to solve geophysical flow problems by using SPH, that groundwater flow, pore pressure changes and precipitation-induced infiltration are major factors causing these geohazards if unchecked.

Keywords:

• Geophysical Flow
• smoothed-particle hydrodynamics (SPH)
• landslides
• debris flow
• geofluid
• slope failure
• mathematical technique
• coupled fluid-solid phase (CFSP)

1. Introduction

It is not an arguable fact that the routine geotechnical engineering and earthwork designer and practitioners do their works is based on the assumptions and field conditions of saturated soils (Houston, 2019). The assumptions that there are no future potential soil wetting, which allows for the use of total stress and the assumptions of total wetting which allows saturated soil parameters to be used for field unsaturated conditions are two major and common approaches to unsaturated soil practice (see, Figure 1; Houston, 2019). It clearly has the potential to be un-conservative to assume no wetting conditions, in which case there is no change in soil suction (Houston, 2019). Conversely, total wetting assumptions are overly conservative, which results in excessive cost implications. Obviously, the difficulties encountered in managing and designing for unsaturated soil conditions, which are time and parameter evaluation and measurement and determination of soil suction can be overcome if soil behaviour problems are to be solved sustainably and efficiently. It is even a complex exercise to design for the safety of slopes under these soil behaviour conditions. Fluid and particles tend to flow under gravity provided they assume height under predefined gradients (Forrester, 1996). This of course is the origin of geophysical flows, which takes the form of landslides, debris flows, slope failures, erosion, etc. (Forrester, 1996). It is important to note also the technical and design difficulty involved in understanding the behaviour of the water table in a sloppy earth arrangement. In this condition, water table becomes a function of the distribution of the particles and infiltration, which relies on permeability and the geometry of the slopes as presented in Figure 1. Geophysical flows (failure and deposition down a slope) happen for various reasons; water in a slope (geofluid), which reduces the stability by reducing the shear strength and increasing pore water pressure (u = z ${\gamma }_w$). This encourages particle dislodgement, deposition, and exertion of water pressure on tension cracks, which results in softening as the water content increases, piping which is caused by groundwater flow (seepage) and changes in the mineralogical composition of soils. Precipitation-induced infiltration and groundwater (geofluid) are the major causes of slope failures due to increased pore pressures (u = z ${\gamma }_w$) and erosion of internal and external slope materials resulting in geophysical flows as shown in Figure 2. Groundwater is considered as stress degradation, which confront slopes and geotechnical structures (Forrester, 1996; Nooralddin, 2021). Slopes are bound to fail after a long period of rainfall as groundwater flows from slopes to some extent (see, Figure 1) depending on the draining rate of the soil, and this causes failure when discharge rate out of the soil mass is not sufficient, which grants excess pore water pressure (u = z ${\gamma }_w$) build-up (González de Vallejo & Ferrer, 2011). In order to evaluate the adverse state conditions of the groundwater, appropriate stability analysis methods are important to estimate the effect of the water flow on the stability of the slopes and on the hydrogeological behaviour of the environment under study. Shear strength of failure surface, (τ = c + (${\sigma }_n$- u) tan φ = c + ${\sigma }_n^{\prime }$ tan φ) reduction due to reduced normal stress as a result of pore pressure build-up (Sharma & Bora, 2003), pore pressure exerted in tension cracks causing down-slope shear forces and soil material weight increase due to saturation (γ = ${\gamma }_d$ + Sn ${\gamma }_w$) and the major hydrogeological considerations and these result to underground flow and softening (Forrester, 1996). The stability analysis of geophysical flows to solve landslide and slope failure problems is a continuous but discontinuity state conditions exist known as tension cracks and an approximate approach to evaluating the force exerted by geofluids on discontinuity surfaces is to assume triangular hydrostatic pressure distribution on the surfaces as illustrated in Figure 3 on the slope and slope head. However, coupled fluid-solid phase state condition is yet valid for stability analysis.

Figure 1. Infiltration and geofluid flow through saturated and unsaturated zones of a slope (Forrester, 1996; González de Vallejo & Ferrer, 2011; Nooralddin, 2021).

Figure 2. Effect of infiltration and soil particle arrangement in water table in a slope (Forrester, 1996; González de Vallejo & Ferrer, 2011; Nooralddin, 2021).

Figure 3. a) Water pressure triangle for a discontinuous single plane on a slope. b) Pressure triangles through tension cracks at the slope head (U and V are forces due to water on the sliding plane and on the tension crack, respectively; Forrester, 1996; González de Vallejo & Ferrer, 2011; Nooralddin, 2021)

Geophysical flows are coupled fluid-solid phase (CFSP) state conditions, which are analysed using complex analytical and/or numerical techniques (Lind et al., 2020; Franci et al., 2020a &, 2020b; Nonoyama et al., 2013). SPH model can accommodate the spatial variability of the CFSP state conditions experienced in landslide, debris flow and slope failure problems due to its mesh-free versatility and powerful flexibility (Lind et al., 2020; Wang et al., 2016). Figure 4 shows the schematic arrangement of the result function and the process of computation known as SPH convolution when applied in the stability analysis of earth and environmental continuum problems (Bui & Nguyen, 2017; Nonoyama et al., 2015). It can be observed that a target particle, “i” is mapped within a domain of multiple particles “j” known as neighbours (Cercos-Pita et al., 2016; Nonoyama et al., 2013). Figure 5 presents the application of a kernel function in reconstructing a continuous function in SPH modelling (Chen & Qiu, 2011; Mocz, 2011). In this research paper, the sufficiency of a more robust field investigation and soil behavior and climate factors profiling have been emphasized on and this is to have a more effective deployment of SPH constitutive modelling technique in solving geophysical flow problems around the world. This approach is aimed at supporting the findings of Stefan Van Baars’ work, which has suggested that friction angle was sufficient to solving geophysical flow problems. In addition to the use of expanded soil profile in the modelling of slope-driven geophysical flow problems, this overview has suggested the application of a particle-based numerical approach, which applies the Lagrangian solver called the smoothed particle hydrodynamics (SPH). This is fundamentally due to the ability of the SPH to deal with mesh distortion problems associated with mesh-based numerical solutions and its ability to execute the modelling operation without the evaluation of the Dirac Delta Function (DDF), which form the basis of the other mesh-free techniques.

Figure 4. Schematic arrangement of an SPH convolution; result function and the process of computation (Chen & Qiu, 2011; Mocz, 2011; Nonoyama et al., 2015).

Figure 5. Application of kernel function (Gaussian-like curves) in reconstructing a continuous function in SPH (Chen & Qiu, 2011; Mocz, 2011; Nonoyama et al., 2015).

2. Background to the study

Landslide is the movement of soil and rock mass down a slope, which denotes mass wasting under the direct influence of gravity. Primarily, gravity is one of the main causes of landslides, acting on weakened materials that make up a sloping area of land. In some cases, landslides occur slowly over time, that is, land movement on the order of a few meters/yards per month. Hence, the most destructive type of landslide happens suddenly after a triggering event such as heavy rainfall or an earthquake. Generally, it is established that physical and anthropogenic factors are the two main causes of landslides (Skilodimou et al., 2018). In practical terms, almost every landslide has multiple causes just as the cases of landslides in the Gutland Southern region of Luxembourg. The region as presented in Figure 6 with its geographical information presented in Figure 7 is characterized by weak Keuper-Rhaetian clay formation and this soil formation is responsible for two-thirds of the landslide occurrences in the country. Based on the previous studies, the clay soil across this region is erratic and requires rigorous research investigation to arrive at a sustainable solution. The conventional method of limit equilibrium proposed by researchers like Bishop, (1955), Janbu, (1954, 1957) Spencer, (1967), whose principles are based on the assumptions regarding the inclination and location of the interslice forces in the slope through the application of factor of safety will not suffix to the slope instability problem in Gutland Southern region. The fact remains that there is no standardised process to measure slopes factor of safety (FOS). Global concepts are randomly used to provide some degree of safety, mainly driven by a professionals’ judgment or previous experience in deterministic geotechnical design. Too high FOS, although conservative, is usually not cost-effective whereas a too low FOS will result in casualty and economic consequences.

Figure 6. Case study area-Luxembourg.

Figure 7. Geographical information of the case study area-Luxembourg.

The clay soil of Gutland Southern is naturally variable because of its particulate characteristic and formation processes. After the soil initial formation, they are modified continuously by external stresses and various processes such as physical weathering, chemical reactions and decomposition, the introduction of new substances, and in some cases, human intervention (e.g., soil improvement techniques, excavation, filling, induced groundwater variation). As such, Gutland Southern soil profiles are inherently heterogeneous causing highly variable properties both horizontally and vertically (Hight & Leroueil, 2003; Uzielli et al., 2006).

Due to the particulate characteristics of Gutland Southern soil as well as other multi-factors, advanced methods of slope stability analysis such as numerical modelling, constitutive models as well as artificial intelligence-based predictive models are supposed to have many advantages over the limit equilibrium method. The conventional numerical modelling based on different methods, such as finite element method (FEM), boundary element method (BEM), and finite difference method (FDM) has become readily acceptable for slope stability analysis, particularly in situations where the failure mechanism is not controlled completely by discrete geological structures. Widely, many researchers have used this technique (Sarkar et al., 2012; Singh et al., 2008; H. Zheng et al., 2005) to simulate the slope stability with reliable results. Also, a vast amount of research work done by Pradhan et al., (2011) pointed out many successes for engineering application of numerical simulation in slope stability analysis. Furthermore, several numerical studies have been carried out on slope stability failures induced by suction, changes in the groundwater table, and rainwater infiltration of slope (Gasmo et al., 2000; Ng et al., 2001).

Considering the geological formation of the Gutland clay soil, it is imperative to conduct and monitor the slope of this region over time. To effectively develop constitutive models as well as artificial intelligence-based predictive models for the slope stability challenges facing the Gutland region, there is a need to conduct field and laboratory tests to determine material properties of this weakened clay using static penetration tests (SPT) and cone penetration tests (CPT), which evaluates shear strength modulus, i.e., angle of internal friction and cohesion of less than 8° and 13 kPa. Furthermore, tensiometer to determine the clay suction over time, sensors and microwave transmitter to soil moisture variation and soil movement. Based on the landslide problems with the associated challenges in Gutland an extensive investigation through computational soil mechanics and constitutive soil modelling is expected to be developed to foster conservative options that could trigger landslide through severity and sensitivity of soil discrete parameters.

3. The Concept of Smoothed-Particle Hydrodynamics (SPH)

3.1. Smoothened Particle Hydrodynamics for Modelling Geophysical Flows

Geological hazards like debris flow, slope and gully failures and landslides are geological hazards associated with high social and economic risk to affected communities. These hazards are characterized by large post failure deformations, high flow velocity (surface and subsurface), long run-out distances and considerable influence zones (Z. Zheng et al., 2020). Geophysical flows like landslides and debris flows cause huge socio-economic losses to the affected communities. Measures to combat these hazards, design protective structures, characterize possible secondary effects like landslide dams or waves, assess risk and plan for their mitigation is mandatory in the regions susceptible to these hazards. A number of models with varying levels of accuracy are used to quantify susceptibility and assess risk of landslides involving long run-out distances and are broadly categorized as empirical models and rational models (Fell et al., 2008). Empirical models use either methods involving volume change (Cannon, 1993; Fannin & Wise, 2001) or the geometrical relationship between the slope and landslides (Corominas, 1996; Crosta et al., 2003; Hungr et al., 2005) and help in estimating the run-out distance but do not provide information on the velocity, direction and thickness of flow. Rational methods employ mathematical models expressed as partial differential equations and provide quantifiable measures of parameters (Braun et al., 2017; Pastor et al., 2014). Numerical modelling can be useful to estimate variables like pore pressure and rheological parameters through back analysis which are difficult to measure through experiments or in the field (Cascini et al., 2014). Rational models can be classified as discrete and continuum models. Discrete models are best suited for simulating avalanches (Pastor et al., 2014) but are not suited for flow as the rheology of flows cannot be adequately described by colliding particles alone (Pastor et al., 2014; Z. Zheng et al., 2020). Coupled models can be used but computational cost is higher. Continuum models are based on continuum mechanics and accommodate the coupling of hydraulic, mechanical and thermo-mechanical behaviour.

Continuum models are further divided as (i) 3D models based on mixture theory (Quecedo et al., 2004; ii) velocity pressure models (iii) depth integrated models and (iv) simplified depth integrated models such as infinite landslide model. 3D models can represent all the phases (solid, liquid and gas) of a flow but are complex and need great effort in developing codes (Pudasaini & Hutter, 2007) and are the models of the future. Velocity pressure models can characterize the triggering mechanism of geophysical flows that can provide the initial conditions and therefore are used in propagation analysis (Pastor et al., 2014). Majority of the geophysical flows likes landslides and debris flows are fast propagating in nature and hence depth integration approximation can be used to model them. Basal pore pressure can also be included in the depth integration model (Iverson & Denlinger, 2001; Pastor et al., 2009) and two-phase models can also be developed (Pelanti et al., 2008). But pressure and force obtained on structures are not accurate owing to limited data on vertical profile of velocity. This limitation can be overcome by coupling 2D model with local 3D model (Pastor et al., 2014). Depth integrated models provide the best compromise between accuracy and computer time and can therefore be extensively used for modelling geophysical flow (McDougall & Hungr, 2004).Geophysical flows associated with free surface, moving interfaces, deformable boundaries, crack propagation and large deformation are difficult to be modelled using grid-based methods as they result in distortions and numerical complications (Huang and Dai 2014). Lagrangian discretization is an attractive alternative to mesh-based methods. Mesh free methods like material point method (MPM) and smoothened particle hydrodynamics (SPH) method can be used to simulate slope movements with large deformation (Dutto et al. 2105, Yerro et al., 2018) with suitable constitutive laws. SPH, a Lagrangian particle-based mesh free kernel-based interpolation technique that can be used to model geophysical flows based and provides advantages like avoiding mesh distortions and less of computer time (Braun et al., 2017; Huang & Dai, 2014). A particle is represented as a finite volume in continuum scale. SPH was introduced by Lucy (1977) and Gingold and Monaghan (1977) for 3D open space astro-physics-related modelling at astronomic scale and was later adopted to model the propagation of landslides and debris flows by McDougall and Hungr (2004), Bonet and H. Zheng et al. (2005) and McDougall and Hungr (2004) and thereafter hydromechanically coupled models of landslides were introduced by Pastor et al. (2009) which was further applied to specific landslide case studies by Cascini et al. (2014) and Pastor et al. (2014). SPH simulates the entire landslide process starting from initial condition of the slope, landslide initiation and transport and deposition of the soil mass and does not require the assumption of small deformation (Li et al., 2019). SPH simulates geomaterials as particles and has been successfully used model to model not only geophysical flows like landslides (Huang et al., 2012; W. Zhang et al., 2019, 20) but also embankments (Bui & Fukagawa, 2013 and river levees (Maeda et al., 2006)as a variety of soil models have been implemented into SPH (McDougall & Hungr, 2004; Pastor et al., 2009; An et al., 2016; Mao et al., 2017; Li et al., 2019; Z. Zheng et al., 2020). Three-dimensional simulation of large deformations has also performed using SPH.

McDougall and Hungr (2004) introduced a simulation model based on depth integrated Lagrangian solution for rapid landslides across terrain modelled in 3D based on the work of Hungr (1995) with depth integration adopted normal to sliding surface. The model uses a meshless technique where mass conservation is satisfied by interpolation and is based on SPH. The volume of the sliding mass is divided into elements called particles and each particle has a finite volume with constant density that increases on entrainment alone but remains at the centre of moving reference columns and the flow depth is proportional to material volume in the area depending on the proximity of particles. The depth of the flow in the reference column “j” can be calculated as the summation of interpolant (Wang & Shen, 1999);

${\sum }_{i=1}^{n}ViWij$

Where i = 1 to n are the particles in the proximity of radius of influence, W—interpolating kernel and V—volume of particle

The model was tested in closed form solution with many controlled laboratory experiments and a full-scale rock avalanche and yielded promising results. Simulation was run through a series of flume experiments for checking the consistency in performance of the model. In the entire analyses carried out by McDougall and Hungr (2004) only two unknown parameters were calibrated—basal friction angle and angle of internal friction and other parameters were kept constant. The calibrated friction angles were close to values obtained from laboratory testing. Their model allowed for anisotropic and non-hydrostatic internal stress, accommodating both material entrainment and option to select various material rheology along the flow path or within the sliding mass. Also, permits large displacements and decoupling or branching of the sliding landslide material. The model also has the advantages like capacity for extensive calibration against landslide events and speedy back analysis. The model used a concept of “equivalent fluid” to represent the heterogeneous landslide material that involves a pair of resistance parameters constrained by calibration and are selected on trial and error based on back analysis of landslides in full-scale. This approach reduces the dependence on laboratory-based material data and the need for a valid constitutive relationship for full-scale landslide events, allowing for comparison of different rheological models and choice of new ones. After calibration with a number of landslide events, the model is capable of predicting run-out. McDougall and Hungr (2004) concluded that though SPH proved as a powerful tool for analyzing landslides dynamically further developments were needed for testing and implementation of the material entrainment algorithm. Investigations of the model on irregular paths as in the shock waves were required to be made and calibration of the model to suit all types of all landslides were to be developed.

Bonet and H. Zheng et al. (2005) applied SPH with a modified version of Langrangian shallow water (SW) equation by extending the SPH model developed by Bonet et al. (2004) with an objective of applying to flows over steep non-uniform terrains as in the case of debris flows and avalanches. The SW assumption is based on the 2D plan view where each particle represents a fluid column of a specified height and the particles move with respect to the topography in the direction tangential to the terrain. The fundamental assumption is that velocity is uniform throughout the height of the column and is parallel to terrain indicating that the spatial variation of height is small and the Lagrangian particles follow in time and the motion of the water column is constrained to move according to the terrain. The storage requirements and computational capacity for the model are minimal as the algorithm solves two components of 3D space only. The modified SPH is robust, stable and allows for a better treatment of the terrain making it more suitable for problems associated with breaking dams, tidal waves, debris flows, avalanches and floods. But they also reported that the major limitation of the model was high implementation cost.

Pastor et al. (2009) proposed a coupled depth -integrated SPH model of velocity pressure type (Biot-Zienkiewicz) for simulating flow-like landslides. The prediction capacity of the model was assessed for problems with available analytical solution and through the simulation of the landslide in Sarno-Quindici area (Southern Italy) in May 1998 which had a volume of 10⁵ m³. The triggered landslide had two source areas and joined the main gully after travelling 1.5 km. The landslide was simulated at the propagation stage and the results of the simulation were compared with the field data. The numerical simulation was performed using a digital terrain model (DTM) of size 5 m x 5 m. Frictional rheological model was adopted for the study and the model accommodated the dissipation of pore pressure in the sliding masses. The pore water pressure varied between 0 and 1, where 1 corresponds to liquefaction and the ratio of the width of the saturated layer to the depth varied from 0.25 to 1. In dry condition the simulated run-out distance was observed to be shorter than the observed run-out distances. Much comparable run-out distances were obtained by taking into account the soil-water mixtures and accounting for coupled behaviour between the fluid and solid phases. The results indicate the good performance of the simulated model and also highlight the importance of the rheological properties of the material selected for the model in estimating the velocity of flow, height of the propagating landslide mass and run-out distance. The simulated model showed that the sliding mass joined at the toe and propagated together before coming to a stop after about 60s. The shape, height and extent of the simulated deposition zone matched the real-time observed deposition zones. The velocity of the sliding masses was also in good agreement. The authors observed that the computational time is lesser than the classical approaches like the finite element methods.

Pastor et al. (2014) used a depth-integrated SPH model to simulate the propagation phase and the model incorporated the dissipation of pore pressure in the sliding mass for problems which had solutions like wet dam break model and real events like Tate’s Cairn debris flow, Tsing Shan debris flow, Cougar hill flow slide and Fei Tsui Road landslide. The pore pressure dissipation was introduced in the model as a separate equation to address flows with high mobility on initiation and stops to move when there is total dissipation of pore pressure. The role of pore pressure in aiding the flow depends on both the consolidation of the soil in flow and also the propagation times. Rheological models like Bingham, Voellmy, frictional and viscous cohesive-frictional models that describe basal friction were used and also simple erosion laws were introduced. Langrangian SPH model was used to discretize the equations including pore pressure which was associated with SPH nodes that were moving. The study also describes the inaccuracies resulting from discretization and also points out that calibration to a large extent depends on the quantity and quality of the data available. Depth-integrated models though use limited material parameters that were obtained through back analysis which is a trial-and-error process but consolidation properties and friction angle were determined from laboratory experiments. The authors report that the initial mass and pore pressure are primary parameter that dictate the results and has to be assumed in general as these parameters are not available. The limitation of the study was the digital terrain models (DTM) resolution which in some cases act as a constraint.

Cascini et al. (2014) proposed a quasi-3D coupled SPH model to simulate real events in Southern Italy. Two mountain basins, Toustolo and S. Francesco, were selected for the numerical analysis. The selected areas are characterized by varying hill slope characteristics and trigger mechanisms. Geo-Flow SPH model (Pastor et al., 2009) was used for the study. The model considers the sliding mass as a mixture of saturated solid skeleton and the variables are the velocity of the soil skeleton (solids) and pore pressure. The adopted model takes into account the bed entrainment along the travel path of the sliding mass and there is a consistent decrease in the elevation of ground with time. Erosion laws of empirical nature can also be incorporated into the Geo-Flow model. The model uses a 3 m x 3 m DTM that reproduces the topography of the selected sites before the landslide events. Landslide inventories and soil thickness maps are used to identify the source areas and initial depth of the landslides. Rheological properties were selected from Pastor et al. (2009) model and the initial pore pressure is assumed as 1 in the source area (i.e) normalized to liquefaction pressure. Automatic time stepping is used to discretize time and time steps shorter than 0.8 s is selected with Runge-Kutte algorithm is adopted for numerical integration. The propagation path and of the sliding mass and its run-out distances of simulated landslide events match the real events indicating the success of the numerical models. The results point out that the rate of bed entrainment and extent of areas susceptible to erosion affect the path taken by the sliding mass, its velocity and deposition height. The selected landslides had extensive dataset available for modeling and also the numerical results yielded satisfactory back analysis for parameters consolidation coefficient, friction angle and growth rate of landslide supporting the choice of frictional approach for debris flow modelling (D’Agostino et al., 2013). The initial position of the water table is consistent with the model results. Limitations of the model are the use of empirical entrainment laws, choice of less detailed DTM and lack of full 3D approach that will help to simulate pore pressure variation in sliding mass effectively.

Cuomo et al. (2016) addressed the lack of research contributions in the space-time evolution of erosion which affects the propagation of debris avalanches and flows and evolution of debris avalanches to flows and vice versa. A quasi-3D coupled SPH model was used to study the role of bed entrainment in the propagation of flows. The study focused on three types of propagation (i) debris avalanche—propagating on open slope (ii) debris avalanche turned to debris flow—channelized flow and (iii) several flow-like landslides converging to the same channel. A case study of flow-like landslide in three different areas in Southern Italy was selected for the study due to availability of data like soil properties, field conditions and landslide evidences at the site. A 3 m x 3 m DTM was used to model the selected terrain and a uniform height of 1.5 m was adopted for the source areas with rheological parameters basal friction, relative pore pressure (equal to 1) and water height (0.25). Consolidation coefficient of 10⁻² m²s⁻¹ was adopted for the study. The high consolidation coefficient indicates the capacity for faster dissipation of pore pressure. The empirical erosion parameter was taken as 0.007 based on best fit from simulation. The results of the numerical analysis reproduce landslide propagation effectively. The results show that eroded thickness varies to a great extent spatially and the maximum thickness are observed at intermediate portion owing to steep slope angles and greater velocity of flow. But the variations in time and erosion rates are similar to that observed in literature in case of debris avalanches. Study also points out that results are site specific and not general.

Braun et al. (2017) modelled the velocity, deposition height and run-out distance of the Baishuihe debris flow to assess the damage and possibility of landslide dam formation using depth integrated hydro-mechanical coupled SPH model using a 5 m x 5 m DTM. The geotechnical parameters were estimated using back analysis due to limited availability of in-situ data. Details on the volume of the landslide, deposition height and its mechanical properties were obtained through field investigation and were used for the back analysis. Step-wise parameter optimization was carried out to reconstruct the observed deposition scenario. The initial volume of the landslide was 0.5 million cubic meters and the triggering height was 15 m. The initial pore pressure was set as 0.6 of liquefaction value with bed entrainment. The height of the water table was considered as half that of soil thickness. The results show that the model error was negligible and the deposition rate was marginally overestimated when compared to the observations in the field. The analysis did not take into account the multiple features involved in the landslide process like initial conditions existing in the source areas, the pattern of flow/propagation, mobilized total volume and features of deposition which could have improved the results of the model making them suitable for assessing the potential of landslide to block rivers.

Li et al. (2019) used the mesh free SPH model to simulate landslides with large displacements after initiation and performed a series of stability analyses on clay slopes using limit equilibrium method (LEM) and SHP which showed that the assumption of using critical slip surface as failure surface is not suitable for unstable slopes. LEM slip surfaces differ significantly from the failure slip surfaces and are therefore not safe for design purposes. The new method of combining LEM with SHP to identify failure surface showed that slip surface with a factor of safety of 1 corresponds to failure in lieu of ones with minimum factor of safety as obtained from conventional LEM.

3.2. Efficient and Sustainable Application of SPH in Solving Geophysical Flow Problems

SPH is a versatile mesh-free numerical simulation technique that is both powerful and flexible for modelling geophysical flows. It has the capacity to predict free surface flows that show fragmentation without loss in mass and is capable of handling large deformations, which are characteristic of geophysical flows like debris flows, rock avalanches, debris avalanches, mudflows and landslides (these are all slope-driven failure). The SPH model can accommodate the spatial variability of the geomaterials and the uncertainty associated with it. It is suitable to model pore pressure variation in the sliding mass, effect of bed entrainment and a variety of rheological features can be adopted to model the geomaterial in the sliding mass. SPH can be applied to varied terrains using suitable digital terrain models (DTM) to model the topography of specific regions though at times the resolution of the DTM presents a constraint. In contrast to the limit equilibrium (LE) and mesh-based numerical methods like the FEM, DEM, FDE, etc., SPH is suggested as the decisive model technique due to its ability to handle distortions that exist in the operations of the discretized meshes of the flow domain of large deformation problems. Also, contrasting the SPH superior abilities with those of the other mesh-free techniques like the MPM and BEM, it’s modelling path is straightforward as a Lagrangian-based application while the MPM and BEM are based on Euclerian biases. Their mathematical formulations evaluate the Dirac Delta Function (DDF) prior to the determination of the final model, which makes the model operation more cumbersome. The numerical approach is reliable, robust and provides an excellent combination of accuracy and computational cost efficiency with good prediction capacity. Generally, for the SPH technique to successfully model the geo-flow problems being studied, the climatic factors other than friction angle like the rainfall intensity, groundwater level (GWL), pore-water pressure (PWP), runoff rate of the watershed, slope changes from the vertical place, intergranular force of the studied soil, frequency of the rain and peak periods, etc. are studied. Additionally, the previous research papers which dealt with the application of SPH included various climatic and geotechnical functions in their analysis of the problem solutions, however, none included the effect of slope angle, heterogeneity of the problem domains and the intergranular forces that exist between particles. It is suggested also that the effect of these missing factors on the failure of slopes in the form of geophysical flows be incorporated in the SPH solution interface.

4. Conclusions

From the foregoing review relevant to the application of smoothed particle hydrodynamics (SPH) in the modelling of geophysical flow problems around the world with emphasis on the central and southern regions of Luxembourg lying within the geographical location presented in the GIS map, which is prone to geohazards, the following can be remarked;

Geohazards are classified by large post failure deformations, high flow velocity, long run-out distances and considerable influence zones and they cause huge socio-economic losses to affected human communities that measures to combat the hazards, design mitigating structures, characterize possible secondary effects, assess risks and plan their mitigations are mandatory in the regions susceptible to these hazards (landslides, slope and gully erosion watershed failures and debris flows). The foregoing has shown that friction angle as a soil variable proposed by Stefan Van Baars can never be sufficiently used in solving landslide problems by adopting evolutionary models like the proposed SPH constitutive relation. It has further shown that for SPH to be effectively deployed, a more robust field study is needed to measure all the soil behaviour and climate factors. A number of models with varying levels of accuracy are used to quantify susceptibility and assess risk of landslides, debris flows and slope and gully erosion watershed failure involving long run-out distances and are broadly categorized as empirical and rational models. Empirical models involve volume changes and geometrical relationships and help in evaluating run-out distances but do not provide information on velocity, direction and thickness of flow while rational models employ mathematical models expressed as partial differential equations (PDEs) and they provide quantifiable measures of parameters. While numerical models can be useful to estimate variables like pore pressure and rheological parameters through back analysis, which are difficult to measure through laboratory and field experiments, rational models can be classified as discrete and continuum. Discrete models are suited for avalanches but are not suited for flows as the rheology of flows cannot be adequately described by colliding particles alone. However, while coupled models can be used in this case, but computational cost is higher. Continuum models are 3D (accommodating solid, liquid and gas states), are based on continuum mechanics and accommodate the coupling of hydraulic (fluid flow under force), mechanical and thermo-mechanical behaviour hence the successful use of SPH in the modelling of geophysical flow problems under coupled situation. Generally, SPH is a mesh-free constitutive technique and does not cause distortions and numerical complications in computation hence its use over other constitutive relations in geophysical flow problems. Therefore, it has the potentials of modelling and solving geophysical flow problems in Luxembourg but sufficient field data and profiling have to be provided beyond what has been suggested by Stefan Van Baars in his work and also this technique can be applied in other parts of Europe, Asia, Africa and America where geophysical flow-susceptible sites exist.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data Availability Statement

The underlying data supporting the results of this study has been reported in the manuscript.

Additional information

Funding

The authors received no direct funding for this research