A Lagrangian–Eulerian procedure for the coupled solution of the Navier–Stokes and shallow water equations for landslide-generated waves

This work presents a partitioned method for landslide-generated wave events. The proposed strategy combines a Lagrangian Navier Stokes multi-ﬂuid solver with an Eulerian method based on the Boussinesq shallow water equations. The Lagrangian solver uses the Particle Finite Element Method to model the landslide runout, its impact against the water body and the consequent wave generation. The results of this fully-resolved analysis are stored at selected interfaces and then used as input for the shallow water solver to model the far-ﬁeld wave propagation. This one-way coupling scheme reduces drastically the computational cost of the analyses while maintaining high accuracy in reproducing the key phenomena of the cascading natural hazard. Several numerical examples are presented to show the accuracy and robustness of the proposed coupling strategy and its applicability to large-scale landslide-generated wave events. The validation of the partitioned method is performed versus available results of other numerical methods, analytical solutions and experimental measures.


Introduction
A landslide-generated wave (LGW) [1], also called landslide-induced tsunami [2], is a major natural hazard that occurs when a landslide impacts a water reservoir and produces large-amplitude waves. LGW events can have devastating effects on the coastal areas of water basins, such as lakes, fjords and artificial reservoirs.
The fjord district of western Norway is one of the zones of the world most affected by this major natural hazard [3,4]. Historical records over the last 400 years show that Norway has experienced at least two major LGW events every century [5]. Only in the first half of the last century, the catastrophic events of Loen (in 1905 and 1936) [6] and Tafjord (in 1934) [7] caused the death of 174 people. The LGW events of Lituya Bay, Alaska, in 1958 [8] and Vajont, Italy, in 1963 [9] are among the most well-known cases of this cascading natural hazard. A wider overview of LGW historical events can be found in [10] and [11]. Furthermore, this situation is made even more critical by the effects of global warming, which is clearly leading to an increment in number and intensity of natural disasters [12][13][14].
Accurate modeling and prediction of LGWs are of key importance to reduce their catastrophic effects. Both experimental and numerical studies have been greatly contributed to the enhancement of the forecasting capabilities against these natural hazards.
Physical models are particularly helpful to identify the key parameters of both the sliding material and the water body and to determine their specific effect on the LGW scenario [15][16][17][18][19][20][21]. A detailed overview of experimental tests applied to LGW events can be found in [1]. Nevertheless, laboratory tests are mainly devoted to determining the near-field wave conditions, while estimations on the far-field waves, which are responsible for major damages to the coastal areas affected by an LGW event, are more difficult to extrapolate.
On the other hand, numerical methods have the potential to predict both near-and farfield waves characteristics. However, the numerical simulation of a LGW is a challenging task. Indeed, the computational method must be able to model the complex constitutive behavior of the landslide material, deal with fluid-solid (or multi-fluid) interaction, and track the largely evolving topology of both landslide and water bodies. Furthermore, the LGW analysis involves different characteristic time and space scales for the near field (landslide-water impact zone) and the far field (wave propagation). Finally, it is required to solve large-scale three-dimensional (3D) geometries for long time durations, and this makes the computational cost of LGW analyses one of the most critical issues.
The numerical models applied to LGWs can be classified into three main groups [11] briefly summarized below.
The first approach consists in using a wave propagation solver, typically based on Shallow Water (SW) equations. In this strategy, the landslide runout and water impact are not resolved but are introduced into the model as an equivalent boundary condition [22,23]. This approach is the simplest one and has the lowest computational cost. However, it assumes strong simplifications on both the landslide motion and momentum transfer and thus it can only give an approximate idea of the global LGW scenario.
In the second strategy, the landslide and water motion equations are solved in a unique coupled model. First applications of this holistic strategy can be found in [24] and [25], where shallow water models were used for both the landslide and the water body. Only recently, more accurate 3D monolithic approaches for LGW problems have been presented, see for example [26][27][28][29][30][31]. Nevertheless, the computational cost of this fully resolved method can be still unaffordable for large-scale events.
The third approach splits the LGW problem into two simulations that interact with each other at their interface. Typically, in these so-called partitioned strategies, a numerical method, here called near-field solver (NFS), computes the landslide runout, impact against the water body and wave formation. A different numerical scheme, here called far-field solver (FFS), predicts the far-field wave propagation [11]. Generally, a weak (or one-way) coupling scheme is adopted, meaning that the NFS is insensitive to the FFS solution. The one-way coupling simplification preserves the computational advantages of this partitioned approach and it still ensures an accurate modeling of the key phenomena of an LGW scenario, such as the landslide runout, the wave generation and the far-field wave propagation.
One of the first applications of this partitioned method for LGWs was presented in [32]. In this work, a simplified 3D model was used for the landslide-water impact and a shallow water model was applied for the far-field wave propagation. In [33], a potential tsunami scenario induced by the collapse of a part of Cumbre Vieja Volcano of La Palma island, Spain, was studied by coupling a 3D compressible Eulerian solver with a Boussinesq model. In [34], the same case study was analyzed using a 3D Volume Of Fluid (VOF) method, as the NFS, and an analogous FFS such as that used in [33]. More recently, Tan et al. [35] coupled a Smoothed Particle Hydrodynamics (SPH) method with a shallow water equations solver was used to reproduce hypothetical LGW scenarios at Es Vedrà, Ibiza, Spain.
In this work, we propose and validate a novel partitioned model for LGWs. In this new strategy, a Lagrangian finite element method, namely the Particle Finite Element Method (PFEM) [36][37][38], is used as the NFS and a standard shallow water Boussinesq model is used as the FFS. Several previous works have shown the accuracy of the PFEM to model landslides [39][40][41], also in cascading events [42][43][44][45]. In this work, we use the PFEM approach that has been successfully applied to LGW scenarios in [46] and in [28,29], where 3D simulations of the Vajont disaster were presented. This work aims at being a proof of concept of this new coupled strategy for real LGW scenarios. For this reason, a deep validation of the method is presented by analyzing the performance and accuracy of the new partitioned technique in targeted tests, using reference solutions obtained with other numerical methods, experimental tests and analytical solutions. In partitioned methods, the momentum transfer between the Navier-Stokes and the Boussinesq models must be accurate in order to obtain a faithful representation of the LGW scenario. Thus, particular attention is devoted here to analyze the effect of the near-field boundary conditions on the far-field propagating wave. Convergence and sensitivity analyses are carried out in order to verify the accuracy and robustness of the proposed method.
The content of the paper is structured as follows. The NFS and FFS are presented in "Near-field and far-field solvers" section. In "One-way coupling" section, the coupling strategy between these two solvers is explained. In "Validation tests" section, three LGW examples are analyzed to validate the proposed partitioned strategy and to show its potential to large-scale LGW events.

Near-field and far-field solvers
In this section, we describe briefly and separately the near-field and far-field solvers that are here used to model LGWs. In section Near-field solver: Lagrangian Navier-Stokes model, we present the Lagrangian multi-fluid Navier-Stokes model that is used as the NFS and, in section Far-field solver: shallow water equations, we introduce the Eulerian shallow water Boussinesq solver employed as FFS.

Near-field solver: Lagrangian Navier-Stokes model
In this work, the Particle Finite Element Method (PFEM) is used to model the landslide runout, its impact against the water body and the wave generation. The motion of both water and landslide materials is obtained by solving the Navier-Stokes equations using an updated Lagrangian description of motion. Following other PFEM approaches [47,48], the mass conservation equation is not solved in the standard divergence-free form but admitting a small compressibility (quasi-incompressible formulation). In this framework, the governing equations read where v is the velocity field, σ is the Cauchy stress tensor, b is the body force per unit of volume, p is the pressure, ρ and κ are the material density and bulk modulus, respectively, t is the time, t is the updated computational domain, and T is the total time duration.
The system (1) is complemented by appropriate boundary conditions that read being n is the outgoing normal vector to fluid boundaries,v the prescribed velocities on the Dirichlet boundary ( v ), andt the tractions acting on the Neumann contour ( t ).
To deal with incompressible materials such as water, the Cauchy stress tensor is split into its deviatoric and volumetric parts as where τ is the deviatoric stress tensor and I is the second-order identity tensor. The deviatoric part of the Cauchy stress can be written in a general form as whereγ is the deviatoric strain rate andμ is the apparent viscosity. We remark that Eq. (4) can represent both water and landslide materials. Water is modeled with a standard Newtonian law and, in this case,μ of Eq. (4) is simply the dynamic viscosity. On the other hand, in the non-Newtonian laws for the landslide material,μ may also depend on the deviatoric strain rate, such as in a Bingham law, and on pressure, such as in a frictional viscoplastic model. Implementation and validation of the mentioned non-Newtonian laws in the PFEM formulation used in this work can be found in [49] and [28].
We also remark that, in this work, the sliding body is always modeled as a mono-phase material, also in underwater conditions. Considering the high impact velocities of the landslide, we can assume that this simplifying hypothesis has an almost negligible effect on the wave generation mechanism. On the other hand, this may affect the final deposition pattern of the sliding material, which, however, is not among the main interests of this article.
In the PFEM, the solution of the governing equations is obtained like in a standard Lagrangian Finite Element Method (FEM). Linear shape functions are used for the nodal unknowns of the problem, namely the fluid velocities and pressure. The formulation is stabilized with a Finite Calculus (FIC) method [50] to avoid the numerical instabilities due to the unfulfillment of the inf -sup condition [51] in the incompressible case. An implicit Fig. 1 a Diagram of the variables for the shallow water model. b Mean velocity. c Velocity at a specific depth two-step strategy is used to solve each time increment. Once convergence is reached, the position of the nodes is updated and the quality of the discretization is checked. If the mesh gets too distorted, a new FE mesh is built maintaining the nodes of the previous one. In the PFEM, remeshing is done with a fast algorithm that combines the Delaunay triangulation [52] and the Alpha-Shape method [53]. A more detailed description of the PFEM approach and the remeshing algorithm can be found in the review work [38].
The Lagrangian description and an efficient remeshing strategy, make the PFEM capable of detecting accurately the highly deforming contours of both landslide and water bodies, which is a key feature for a correct simulation of a LGW event.

Far-field solver: shallow water equations
The Shallow Water (SW) equations are a simplified form of the Navier-Stokes equations. The main assumptions include incompressibility, hydrostatic pressure distribution, small vertical velocity. After integrating the equations from the seabed to the free surface η, a new set of unknowns is obtained. The pressure p is replaced by the free surface elevation η, and the velocity u is substituted by the mean horizontal velocityū. We will use the modified Boussinesq equations presented in [54]. Those equations are an extension of the SW equations which include the modeling of the frequency dispersion for long waves. Originally, these equations were expressed in terms of the horizontal velocity u β , evaluated at a specific relative depth β, namely, The auxiliary fields J η and J u introduce the dispersive mechanism and are defined according to the following expressions: where the C i constants depend on the choice of β The free parameter β was fixed to −0.531 by Nwogu et al. [54] in order to minimize the errors introduced by the approximation with respect to linear wave theory. Fig. 1 shows a schematic of the model. System 5 is closed with the following Dirichlet boundary conditions at the reflecting R and inflow I boundaries, Where u and η are the specified inflow velocity and the wave amplitude. Both variables are correlated by linear theory ( [54,55]). In the case of supercritical regime, both velocity and wave amplitude should be imposed. Due to the dispersion relation [56], the reflecting boundary condition has the following additional requirement if the velocity is imposed The Boussinesq equations are solved using the standard FEM. The space domain is interpolated with a Galerkin discretization of linear triangles and a finite difference scheme with constant time step is used to integrate the equations in time. Some numerical difficulties such as the third-order differential operator and the time integration accuracy are addressed in [56,57]. As stated in [58,59], problem (5) is an hyperbolic wave in mixed form and there is an incompatibility condition (see [51]) because the same interpolation is used for both variables, the velocity u β and the wave amplitude η. Here, the equations are stabilized using the FIC approach extending the work reported in [60].
Following [57], the third-order spatial derivatives are modeled using J η as an intermediate variable and the field J u can be directly included in the equations by parts integration. Some boundary terms arise from the integration of both fields. Since these terms do not vanish at all the boundaries, we cannot neglect them.
The time integration has been approximated using the predictor-corrector Adams Bashforth Moulton method (see for example [56,61]).

One-way coupling
The LGW problem is here simulated using a weakly coupled (one-way) method which makes interact the PFEM solver presented in section Near-field solver: Lagrangian Navier-Stokes model, and the SW solver described in section Far-field solver: shallow water equations.
We note that problem (5) defines a phase speed c which in the SW limit is computed as gh. Being u the modulus of the horizontal velocity u β , the information travels at velocities u + c and u − c. If the flow is subcritical, which is the case analyzed in this work, then u < c and the information travels both upstream and downstream.
Considering the bidirectional characteristics of the equations, a first possibility would be to consider two domains adjacent to each other and define a strong (two-way) coupling (see for example [62]). However, such a strongly coupled approach is computationally expensive, as it requires running parallelly the PFEM and the SW solvers.
Moreover, the accuracy increment given by a two-way method over a one-way strategy on the far-field wave propagation can be considered negligible. For these reasons, here, we adopt a decoupled space-time one-way strategy in which the PFEM solution is stored at the SW interface and, in a second stage, it is imposed as a boundary condition for the SW simulation. In order to avoid perturbations of the results at the interface, the computational domain of the PFEM is extended beyond the position of the SW interface using non-reflecting boundaries. Taking advantage of such a one-way coupling strategy, the PFEM and SW simulations can be executed independently leading to a very versatile tool for LGWs with significant saves of computing time. Since the PFEM is a Lagrangian strategy, a search algorithm is constructed at every time step in order to find all the elements cut by the SW interface. Then, the PFEM calculations beyond the interface are not relevant. This fact is the key to the computational savings, since the computational domain can be shortened by means of an open boundary. However, the numerical approximation of open boundaries-the absorbing boundaries-introduces some reflections. In this work, the absorbing boundary is modelled by extending the domain after the open boundary with a gentle slope. The computational domain ends when the slope reaches the mean water level, at this point, the impulse waves leave the computational domain.
In a later stage, the characteristic variables computed at the interface are imposed to the SW domain through an inflow boundary condition. We recall the subcritical characteristics of the analyzed flows, hence, one variable is required to be imposed in order to define a well-posed problem: the wave amplitude or the horizontal velocity. We choose to impose the velocity, since it is more representative of the momentum exchange from the PFEM and the SW computation. It has proven to be accurate, even when the Boussinesq assumptions are not perfectly fulfilled. A general picture of the coupling strategy is illustrated in Fig. 2.
Even though the Boussinesq equations are expressed in terms of the velocity evaluated at a certain depth, u β , this magnitude is a measure of the depth-averaged velocityū. In other words, it can be understood as a numerical quadrature of one integration point. When the waves are regular, the choice of one magnitude or another is not relevant, but when wave breaking is present, the depth-averaged velocity is more representative of the momentum exchange.
We remark that the average vertical velocity of the fluid corresponds to the time derivative of the free surface elevation. This variable does not correspond to a boundary condition for the studied cases.
Finally, there is an additional condition associated to I (8b). Ifū is assumed to be equal to u β , a constraint appears over the dispersive field J η , which is equivalent to impose ∇∇ · u = 0 on I (see [57]).

Validation tests
In this section, we present three different cases to validate the proposed partitioned strategy and to show its potential for practical applications. The first numerical example is aimed at reproducing a unidirectional wave generated in a laboratory channel. For this test, we carry out a detailed validation of the coupled method paying special attention to the transmission of boundary conditions between the near-and the far-field solvers. The simplicity of this test allows us to compare our results with both experimental measures and analytical solutions, and also with the numerical solution obtained with a full PFEM model. In the second example, we apply the partitioned method to a more representative example of LGW problems. In this test, we reproduce numerically the water wave generated experimentally by the impact of a second mass of water sliding at high velocity over a steep slope. The last test aims at showing the applicability of the method to real-world LGW problems. For this purpose, we considered a realistic configuration of a LGW event occurring in an alpine lake. Our numerical solution is compared to another LGW solver presented in the literature.

Solitary wave in a channel
In this test, we reproduce the laboratory experiment carried out at a large wave flume of the Coastal Research Center in Hannover. A solitary wave is generated by a pistontype maker and travels 180m until reaching the final inclined slope. A schematic view of the wave flume is depicted in Fig. 3. More details about the experiment can be found in [63][64][65]. Figure 4 shows the horizontal stroke of the paddle along time. The wave height has been monitored at different positions of the flume, including the on-shore zone. In this work, we will compare our numerical solution to the experimental measures obtained at the four wave gauges whose coordinates are given in Table 1. The selected gauges are placed   [63][64][65] at key positions of the channel and allow us to monitor wave generation (G1), propagation (G2), shoaling (G3) and flooding (G4).

Physical considerations
The aforementioned specifications generate a solitary wave of 0.6 m amplitude and 65 m wavelength. The wave generation, propagation and breaking were analyzed using the PFEM approach reported by Oñate et al. [66]. Given the properties of such a solitary wave, it can be simulated using the Boussinesq approximation and thus reducing drastically the computational demand. This experiment is very interesting for two reasons. Firstly, we can perform a verification test of both formulations and compare the numerical results against experimental data. Secondly, the simplicity of the geometry allows us to obtain analytical solutions for the Boussinesq equations. The analytical solution is a wave equation of the type Details of the parameters A 0 , A 1 and A 2 and the relation between the wavelength, period and amplitude can be found in [61]. The generation of solitary waves has motivated several discussions and a review can be found in [67]. The kinematic description of the piston wave maker is the origin of the discussion, since it cannot represent the exact solution of a solitary wave due to construction limitations. Some expressions for the motion of the piston can be obtained by integrating the analytical solution of the wave and truncating it on a finite space and time domain, corresponding to the features of the piston. Then, the experimental solitary wave is generated with a tail of secondary oscillations.
The Lagrangian formulation of PFEM perfectly tracks the movement of the paddle and thus the numerical simulation reproduces the experimental results with high fidelity. On the other hand, since the Boussinesq equations are implemented in an Eulerian frame, this boundary condition is difficult to impose. An easier alternative is to apply the analytical solution as a boundary condition.  The analytical solution overestimates the phase speed and this mismatch will be discussed in the following analyses.
We remark that the difference in the phase speed of the wave is not originated by the coupling strategy, but by the Boussinesq approximation. The accuracy of the approximation depends on the non linearity ratio = η/H and dispersion ratio μ = H/λ. A more detailed study can be found in [68], particularly when < 0.4.

Numerical results of the coupled strategy
A global representation of the wave propagation is found in Fig. 6. In this simulation, the first 10 m are simulated using the 2D PFEM and the rest of the channel is simulated using the SW solver. Additionally, the full channel has been simulated with the PFEM to provide a reference solution for the coupled method and to analyze better its performance. Concerning the space and time discretizations used in the two solvers, the PFEM domain has a mesh of mean size x = 0.3 m and the time step increment t = 0.001s is used, while the SW domain is discretized with x=0.8 m and t=0.025 s.
The results of gauges G2 and G3 show a small gap between the predicted wave by the two solvers. The Boussinesq approximation is triggering this gap, originated by an overestimation of the phase speed. This difference is consistent with wave theory and the current wave specifications. Note that the same gap can be observed in Fig. 5. The run-up (G4) is out of the SW theory assumptions, but still relevant results are obtained.
The magnitude of the computational time saving of the coupled method versus the full PFEM solution is about 95%. These savings will be analyzed in more detail in the next paragraphs. The savings depend on the spatial and temporal domains chosen for the NFS, that have to be carefully designed in order not to introduce additional errors.

Sensitivity to the interface position
The FFS sensitivity has been tested with some SW interface positions, namely at x 1 = 10, 20, 30 and 40 m. One would expect to obtain a more accurate response as the interface is placed further away from the paddle. Nevertheless, since in this example the wave is very  regular, the observed influence of the interface position on the results is not significant. The solutions obtained with all the interfaces can be considered already converged (Table  2). These results were expected due to the regularity of the wave. For this reason, a similar study is also performed in section Wave generated by a water landslide, where the SW interfaces are placed into a more chaotic fluid flow.

Sensitivity to the temporal domain
Part of the saving in computational time comes from reducing the duration of the PFEM simulation up to the minimum time needed. Once the initial impulse has generated the wave and it has been transferred to the SW domain, the PFEM computations do not provide relevant information. From that time on, the initial boundary condition, which corresponds to water at rest, is imposed at the SW domain. This transition in the BC has to be carefully treated in order to avoid unphysical oscillations. A good duration for the transition is half of the period of the current wave.
In this test, we evaluate the effect of feeding the FFS with NFS solutions limited in time. In particular, we considered four PFEM analyses of duration 10, 20, 30, and 40s. Figure 7 shows the time evolution of the wave amplitude at the first gauge. In the graph, we also added dots representing the time instant when one analysis starts to diverge from the rest. It is clearly observed that the four solutions have an identical behavior in the first part of the graph. In particular, even with just 10s of the PFEM simulation, the main wave is well reproduced. Beyond this time, the curves diverge progressively. As expected, a time interval of around 10 seconds separates the consecutive diverging points.

Sensitivity to PFEM domain length
Besides the reduction of the time duration of the analyses, the optimization of the size of the PFEM computing domain can drastically reduce the computational cost of the simulations without affecting the accuracy of the results. For this reason, we analyze here the effect of considering partial PFEM domains of 10, 20 and 30 m length plus an extension acting as absorbing boundary condition, as shown in Fig. 8. The study is carried out for both 2D and 3D PFEM domains. The errors introduced by the effect of shortening the PFEM domain are listed in Table 3.
It is important to note that the vicinity of the absorbing boundary condition of the PFEM may affect the accuracy of the interface. The small errors obtained when the interface is far enough from the absorbing boundary show that the presented methodology allows to effectively reduce the PFEM domain without virtually affecting the quality of the solution. This is particularly noticeable in the 2D case.
The 3D case presents a similar behavior, but higher errors are observed in the 30 m domain length. However, these errors are more attributable to the capabilities of the NFS for reproducing the fluid-solid interaction at the lateral walls (see [66] for more details) Fig. 9 Landslide wave problem. Setup of the LGW flume for the experimental and numerical analyses than to the coupling strategy. A finer discretization in the PFEM mesh would reduce this bias.

Wave generated by a water landslide
In the second example, we simulate the experiment carried out at the Queen's University landslide flume presented in [69]. In this laboratory test, a mass of water is released from an elevated reservoir and, after flowing downhill over a 30 • slope, it impacts at high velocity the water at rest placed on a 33.8 m-long channel. In the reference work [69], 41 experiments were presented covering a wide range of source volumes and reservoir depths. In [46], a comparison of experimental and numerical results obtained for three different water depths in the channel is presented. In this research, we select the largest volume case (0.45 m 3 ) and water depth (0.60 m). Figure 9 shows the geometry of the experimental setup considered in this work.
We remark that considering a water landslide does not affect the relevance of the test in the field of LGWs. In fact, the phenomena produced by the water runout and impact are totally representative of a realistic LGW scenario with a fast mobilized material. Furthermore, the use of water as sliding material removes the uncertainty related to the rheological properties of the slide and allows repeatability of the test.
The PFEM is used to simulate the water runout, the impact against the water at rest and the consequent wave formation (Fig. 10). Remarkably, the front of the water landslide reaches the end of the slope with a thin layer of less than 10 cm and it impacts the water in the channel at a speed of about 8.5 m/s. Thus, in order to capture accurately the phenomena at the impact zone, fine mesh and time discretizations are necessary. For this reason, a mesh size of x=1.5 cm and a time step increment of t = 5 · 10 −4 s are used in the PFEM simulations. On the other hand, much coarser mesh and time discretizations can be used to model the wave propagation along the channel with the SW solver. In particular, in the FFS a time step of t = 0.025s and a mesh size of x=0.3 m have been used. We remark that the possibility of using much different and yet adequate space and time parameters in the FFS and NFS solvers is one of the main advantages of this partitioned method and one of the reasons for its high computational efficiency.

Numerical results
This LGW scenario has been solved using a time and space reduced PFEM domain in combination with three SW interfaces. The PFEM spatial domain includes the runout, the first 7 m of the flume and an absorbing boundary condition, while the temporal domain includes only the first 5s. The SW interfaces are positioned at 2, 4 and 6 m. Figure 11 presents the results obtained at the gauges and a representation of the wave propagation. In the top image, it can be observed the vicinity of the first gauge and the first SW interface to the wave generation zone. Indeed, gauge G1 can only record the PFEM solution and the SW solution obtained by the first interface. It is clear that the imposed boundary condition does not satisfy the Boussinesq assumptions and the interpolated wave does not fit the profile of a breaking wave. However, although the wave interpolated by the FFS at the first stages is not equivalent in terms of wave height, the stored momentum is the correct one. This can be observed at gauges G2 and G3, where the experimental wave has adopted the solution of a solitary wave and matches the profile of the FFS. The results obtained at gauges G2 and G3, placed at the middle and the end of the channel, respectively, show that all the three SW interface positions reproduce well the main wave obtained experimentally. This is particularly remarkable considering that the SW interface placed at x=2 m is completely inside the impact zone (Fig. 10). These results show that, as long the momentum is well transferred from the NFS to FFS, the wave propagation process in the far-field can be accurately reproduced even considering the SW interface in a zone where the wave is not completely generated. We also remark that this can be done safely in this test, since water has been considered for the sliding material. In case of considering a different landslide material, either the interface is placed further the zone of material deposition of the landslide, or the interface boundary conditions have to take into account the presence of different materials in the computation of the overall momentum.
Gauge G2 also records a considerable time interval after the first wave, this allows us to analyze also the secondary waves. In this case, we note some discrepancies between the results obtained by three SW interface positions. In particular, the first solution that diverges from the experimental one (and from the two other numerical solutions) is that obtained by the farthest interface position (6 m). This result is totally consistent with the time domain truncation explained in section Sensitivity to the temporal domain and Fig.  7. As the interface position is further from the impact zone, the signal arrives later. Being the phase speed about 2.5m/s, the time difference between each interface is around 0.8s.
As a concluding remark for this example, the computational cost of the full simulation of the LGW has been estimated proportionally to the time needed by the signal to arrive at the end of the channel and proportionally to the number of elements required to discretize the full domain. The resources consumed by the FFS can be neglected since they are two orders of magnitude smaller. According to these considerations, the overall time saving given by the proposed partitioned strategy is 95%. In [70], different metrics of real alpine lakes were used to define the configuration of theoretical mountain basins of different sizes and shapes. These geometries were used in [70] to study LGW scenarios with a finite volume solver and to obtain correlations between the lake configuration and the landslide-generated waves. Here, we analyze one of the lakes considered in [70] to test the proposed coupled strategy in a 3D complex setup. Figure 12 shows the side and top views of the geometry of the lake. The case study is a circular lake with a diameter of 1500 m. The landslide has a prismatic shape of 20 m thick, 208 m long and 120 m wide. Following [70], a bulk material density of 1620 kg/m 3 is used for the landslide material and an initial velocity of 20 m/s has been prescribed to the sliding body.
Preliminary NFS analyses of the LGW scenario showed that the landslide material reaches a deposition distance of around 350 m. This information is useful to place the SW interface at a position that is not trespassed by the sliding material. For this reason, the interface of the FFS has been placed at 400 m from the center of coordinates, which is the center of the run-out impact. Figure 13 shows a global view of the simulated LGW and a superposition of the NFS and FFS results.

Numerical results
In order to assess the quality of the obtained solution, in Fig. 14 we compare the envelope of the maximum wave height measured along sections S1 and S2 with the reference solution given in [70]. Globally, the results obtained with the proposed method agree well with the reference numerical solution, both in the near and far fields. Although with some differences in terms of magnitude, both methods are also able to reproduce the amplification of the wave near the shoreline. This phenomenon is produced by the combined effect of shoaling and the wave reflection given by the steep bottom surface.
We also highlight that the results of the FFS are in good agreement with wave propagation theory. In an unconstrained plane, the wave amplitude is inversely proportional to the distance from the origin. Section S1 is closer to the unconstrained decay, while section S2 shows a smaller decay since it is closer to the boundary.
Finally, it is worth commenting on the peaks in amplitude exhibited by the FFS solution close to the SW interface. As mentioned before, the imposed signal coming from the PFEM simulation is still not fulfilling the Boussinesq theory. On the other hand, the generation of stable waves by Dirichlet boundary conditions requires some traveling distance to be modulated by the fluid system [68]. For this reason, the wave amplitude results obtained close to the SW interface with the FFS should be disregarded. We emphasize again that, on the other hand, the overall momentum computed in that zone is still correct.
In any case, the presented partitioned approach would be really interesting for an exercise like the lakes classification in [70]. Indeed, a single landslide calculated with the NFS could be used to simulate different representative lakes with the FFS. Also, in a more detailed study it would allow concentrating the computational resources in the analysis of the run-out and wave generation, thus enhancing the overall accuracy of the partitioned scheme.

Concluding remarks
We presented a novel partitioned strategy for solving landslide-generated wave (LGW) problems. The coupled method makes interact a near-field solver (NFS) with a far-field one (FFS). The NFS reproduces the landslide runout and the impact zone by solving the Navier Stokes equations with the Lagrangian Particle Finite Element Method (PFEM). On the other hand, the FFS uses as input the NFS results stored at a certain interface to model the wave propagation with an Eulerian Finite Element Method (FEM) based on Boussinesq shallow water (SW) equations. To improve substantially the computational performance of the method and, thus, to allow for the simulation of large-scale problems, we adopt a one-way coupling scheme, meaning that the NFS solution is insensitive to the FFS one. This partitioned method also allows us to freely decouple the time and space discretizations of the two solvers, giving a further advantage in terms of accuracy and efficiency.
In all the examples presented, the results obtained with the new partitioned method had shown a very good agreement with the reference solutions, both in 2D and 3D problems. Remarkably, we have been able to compare our numerical results with analytical solutions, fully-resolved numerical simulations of LGW events, other coupled methods presented in the literature, and experimental observations.
Placing the SW interface as close as possible to the impact zone gives the major advantage of reducing the NFS domain and, consequently, the overall computational cost of the analysis. For this reason, we have compared the FFS results obtained considering different positions of the SW interface for the same NFS solution. This study showed that, as long the momentum of the NFS is well transferred to the FFS, the SW interface can be also placed very close to the impact zone, even if the wave is not already formed. More specifically, the SW interface can be placed at one wavelength from the impact zone. In fact, although locally the FFS results may give spurious amplitudes since the input wave is not fulfilling the Boussinesq theory, the stored momentum is correct and the far-field wave propagation is reproduced accurately. We remark that this can be easily done in case of having the same density between the sliding material and the water in the reservoir, such in the water landslide scenario analyzed in "Wave generated by a water landslide" section. In a more general case, the interface should be placed further than the deposition zone of the landslide, or the SW interface should take into account the variation of material densities on depth.
We have also verified the effect of reducing the size of the PFEM domain by using absorbing boundary conditions. For this purpose, a gentle final slope with an inclination of 1:10 was placed at the end of the PFEM domain. We showed that, as long as the SW interface is not placed too close to the absorbing boundary, the PFEM domain can be safely truncated without affecting the global results. To be precise, the gentle slope should begin at least one half wavelength after the SW interface.
Finally, we also studied the effect of reducing the time duration of the NFS analyses. We have shown that, if the main interest of the simulation of the LGW scenario is to reproduce the main wave propagation, the PFEM analysis can be safely stopped after it has modeled the impact of the landslide on the water and the first wave formation. Indeed, this time truncation of the NFS will only affect the secondary waves propagation. We also showed that, knowing the NFS duration and the wave propagation speed, it is possible to have a quite accurate estimation of the reliability of the secondary waves results.
All these specific studies will allow us the define the most computational efficient NFS-FFS scheme for practical LGW simulations. Although the overall computational cost depends inevitably on the geometry and the proportions among the near and far fields, in the examples here presented, we could estimate a 90% of time saving versus a fully-resolved simulation of the same LGW scenario.
Among the possible enhancements of the proposed method, we consider it of primary interest to investigate more efficient strategies for the NFS absorbing boundaries and to develop a reverse one-way coupled algorithm where the FFS transfers the information to the NFS. This FFS-NFS model would allow us simulating with high accuracy the effect of tsunami waves produced by landslides (or by some other source, i.e., an earthquake) on the shoreline and the civil constructions placed therein.

Funding
The first author gratefully acknowledges the Universitat Politècnica de Catalunya and Banco Santander for the financial support for his pre-doctoral Grant (53 FPI-UPC 2018). This research was partially funded by the project PARAFLUIDS (PID2019-104528RB-I00). The authors also acknowledge the Severo Ochoa programme through the Grant CEX2018-000797-S funded by MCIN/AEI/10.13039/501100011033.

Availability of data and materials
Not applicable.

Declarations
Competing interests