A multilayer shallow water model for tsunamis and coastal forest interaction

ESAIM: M2AN 59 (2025) 1113–1144 ESAIM: Mathematical Modelling and Numerical Analysis https://doi.org/10.1051/m2an/2025017 www.esaim-m2an.org A MULTILAYER SHALLOW WATER MODEL FOR TSUNAMIS AND COASTAL FOREST INTERACTION Raimund B¨ urger1, Enrique D. Fern´ andez-Nieto2and Jorge Moya1,2,* Abstract. Models and numerical methods of the impact of tsunamis on coastal forests are of vital importance for exploring the potential of coastal vegetation as a means of mitigation. Such a model is formulated as a multilayer shallow water system based on a free-surface formulation of the Euler equations for an ideal fluid. Specifically, the Euler equations are approximated by a layer averaged non-hydrostatic (LDNH) approach involving linear pressures and piecewise constant velocities. Furthermore, based on [Iimura and Tanaka, Ocean Eng. 54 (2012) 223–232] drag forces, inertia forces, and porosity are added to model the interaction with the forest. These ingredients are specified in a layer-wise manner. Thus, the vertical features of the forest are described with higher accuracy than within a single-layer approach. Projection methods for the non-hydrostatic pressure in conjunction with polynomial viscosity matrix finite volume methods [Castro and Fern´andezNieto, SIAM J. Sci. Comput. 34 (2012) A2173–A2196] are employed for the numerical solution of the multilayer model, that is for the propagation of tsunamis and coastal flooding. Experimental observations and field data are used to validate the model. In general good agreement is obtained. Mathematics Subject Classification. 65M06, 76D99. Received May 7, 2024. Accepted March 7, 2025. 1. Introduction 1.1. Scope A number of geophysical applications such as shallow water flows, free surface flows, gravity currents, sediment transport and avalanches give rise to a system of first-order partial differential equations (PDEs) in the vectorial form 𝜕𝑡𝑊+𝜕𝑥𝐹(𝑊) + 𝐵(𝑊)𝜕𝑥𝑊=𝑆(𝑊)𝜎′(𝑥),(1.1) where 𝑡is time, 𝑥is the spatial coordinate, and the sought quantity is a vector 𝑊=𝑊(𝑥, 𝑡) of state variables, where 𝑊belongs to an open convex subset 𝐷⊆R𝒩. The vector functions 𝐹:𝐷→R𝒩and 𝑆:𝐷→R𝒩as well as the matrix function 𝐵:𝐷→R𝒩 ×𝒩 and the real scalar function 𝜎=𝜎(𝑥) are given. Solutions of equations of this type are in general discontinuous, and the well-known salient property of Keywords and phrases. Finite volume method, layer averaged non-hydrostatic approach, multilayer model, coastal forest, tsunami mitigation. 1CI2MA and Departamento de Ingenier´ıa Matem´atica, Facultad de Ciencias F´ısicas y Matem´aticas, Universidad de Concepci´on, Casilla 160-C, Concepci´on, Chile. 2Departamento de Matem´atica Aplicada I, ETS Arquitectura, Universidad de Sevilla, Avda. Reina Mercedes No. 2, 41012 Sevilla, Spain. *Corresponding author: [email protected];[email protected] c ○The authors. Published by EDP Sciences, SMAI 2025 This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 1114 RAIMUND B¨ URGER ET AL. the system (1.1) that complicates its analytical and numerical treatment is the presence of nonconservative products such as 𝐵(𝑊)𝜕𝑥𝑊[1,2]. It is the purpose of this contribution to study a specific model of the form (1.1) that arises as a multilayer shallow water system based on a free-surface formulation of the Euler equations for an ideal fluid. The novelty consists in the formulation of this model and its application to forest hydrodynamics as well as its numerical solution by a projection method. Specifically, the Euler equations are approximated by a layer-averaged non-hydrostatic (LDNH) approach involving linear pressures, piecewise constant horizontal velocities and piecewise linear vertical velocities. Furthermore, specific ingredients such as drag forces, inertia forces, and porosity are adopted from the literature of ocean engineering, landscape ecology and related fields [3–7] and are incorporated to form a new model of the interaction of a tsunami wave with a coastal forest. In fact, such models and numerical methods are of vital importance for exploring the potential of coastal vegetation as a means of mitigation. A particular feature of the present approach is the description of tree-specific ingredients in a layer-wise manner. Thus, the vertical features of the forest are described with higher accuracy than within a single-layer approach. The second novel contribution, besides the formulation of the multilayer model, is the formulation of a method for its numerical solution inspired by Chorin’s classical algorithm for the compressible two-dimensional (2D) Navier–Stokes equations [8–10]. In the latter context that algorithm is based on a Helmholtz decomposition of the sought velocity field and proceeds, in each time step, by first updating the velocity field (by using the momentum equation without the pressure term) to give an intermediate velocity, and then by computing the finally updated velocity field by a projection onto the space of divergence-free vector fields. The second step involves the solution of an elliptic problem for the sought pressure. This well-known procedure is followed for the present model, although the formulation of the governing equations is more involved due to the Saint-Venant single-layer or multilayer approach and the presence of forces. Projection methods for the non-hydrostatic pressure in conjunction with polynomial viscosity matrix finite volume methods [11] are employed for the numerical solution of the multilayer model, that is for the propagation of tsunamis and coastal flooding. Experimental observations and field data are used to validate the model. 1.2. Related work Finite volume (FV) schemes are a standard method for solving hyperbolic systems of partial differential equations (PDEs). However, hyperbolic systems arising from balance equations in geophysical applications usually involve nonconservative products that complicate the application of traditional FV schemes. The standard example are shallow water equations with variable bottom topography. A well-known class of FV schemes that can handle nonconservative products are the so-called path-conservative schemes [12– 14]. These schemes are designed specifically for nonconservative systems and are based on the concept of path integration. Path-conservative schemes have been applied to a variety of hyperbolic systems with nonconservative products, including two-layer and layer-averaged shallow water equations, compressible gas dynamics and magnetohydrodynamics. In [15] a generalization of the Roe method [16] was proposed (see also [14]). Nevertheless, its implementation requires explicit knowledge of the eigenstructure of the intermediate matrices. In [11] a specific family of path-conservative schemes is proposed, named polynomial viscosity matrix methods (PVM methods) that generalized several incomplete Riemann solvers such as Rusanov, Lax-Friedrichs or HLL, among others. Multilayer models are designed to avoid solving a fully three-dimensional model (such as the Navier– Stokes equations for an incompressible fluid). They are based on the so-called shallow water or Saint-Venant approach, that is, a vertically integrated version of the underlying model [17–26], in our case the Euler equations for an ideal fluid. The multilayer approach consists in subdividing the computational domain into 𝑁layers in the vertical direction, which leads to a system of Saint-Venant equations. The unknowns in the present case are horizontal velocities by layer, the total height of the fluid column, and pressure. The present approach is based on the LDNH0model with non-hydrostatic pressure. This approach is an improvement compared with standard shallow-water-based models since the latter usually only cover MULTILAYER SHALLOW WATER MODEL FOR TSUNAMIS AND COASTAL FOREST INTERACTION 1115 hydrostatic pressures and neglect vertical acceleration (and therefore dispersive) effects. On the other hand, dispersive models like Boussinesq introduce high-order derivatives for the unknowns, while the LDNH0model incorporates these dispersive effects into the non-hydrostatic pressure terms [27,28]. From the applicative point of view, the general significance of tsumami disaster mitigation by the natural method of coastal forest plantation is discussed in the overview article by Tanaka [6] (for instance). The present work is based on detailed technical information provided in [3–5,7,29]. 1.3. Outline of the paper The remainder of the paper is organized as follows. In Section 2we introduce preliminaries, starting with the basic method for the discretization of (1.1). In Sections 2.2 and 2.3 we study properties the numerical methods need to have in order to be well balanced and behave correctly in dry front situations. To put the LDNH0approach into the proper perspective we first formulate, in Section 2.4, a hydrostatic reconstruction of a two-equation shallow water model with a source term, and then, based on these results, we proceed in Section 2.5 to formulate the linearized discontinuous non-hydrostatic model (LDNH0model) whose unknowns are the height of the water level, the horizontal and the vertical velocity, and pressure as functions of position 𝑥and time 𝑡. The first three of these quantities are specified by a first-order system of balance laws that in each time step can be handled by the same discretization as the shallow water equations along with hydrostatic reconstruction. The non-hydrostatic pressure, in turn, is updated via a projection step that will eventually lead to the final update of all sought variables. (This procedure is reminiscent of Chorin’s method for the 2D incompressible Navier–Stokes equations, as mentioned in Sect. 1.1.) Finally, we specify in Section 2.6 the inertia and drag forces related to the coastal forest along with the corresponding concept of porosity. The result is the non-hydrostatic LDNH0model, specified for forest forces, in final form. Section 3is devoted to the development of the discretization of the LDNH0model. To this end we formulate first, in Section 3.1, a preliminary discretization of the first step of Section 2.5 that excludes the non-hydrostatic pressure terms which are handled by the projection method. This preliminary discretization is based on plausible arguments but turned out to be unstable due to the nature of the inertia force. An alternative, slightly different but stable discretization is advanced in Section 3.2. The numerical method relies on knowledge of the eigenvalues of a matrix related to the Jacobian matrix of the system. These eigenvalues are obtained in Section 3.3. Finally, in Section 3.4 the scheme that discretizes the elliptic problem for the pressure update is formulated. The treatment of Section 3corresponds to the fully discrete version of the projection method akin to Chorin’s method and refers to the single-layer LDNH0model. The multilayer version of that model and its discretization are described in Section 4, starting with the definition of multiple layers (Sect. 4.1). We then present (in Sect. 4.2) the LDNH0multilayer model (without forces and porosity) arising from layer-wise vertical integration. After discussing layer-wise porosity (in Sect. 4.3) we derive in Section 4.4 explicit expressions of the interlayer transfer terms. We then formulate ingredients of the multilayer model that are specific to the application to a forest, namely drag and inertia forces, friction, and viscosity (Sects. 4.5 and 4.6). The resulting multilayer model is summarized in Section 4.7. Next, we outline the discretization of the multilayer model. Roughly speaking, the discretization of the first-order system, described in Section 4.8, is a multilayer version of the discretization of Section 3.2 for the single-layer case (in both cases, non-hydrostatic pressure is disregarded). The remaining ingredients of the multilayer scheme, namely the projection matrix describing the solution of the elliptic problem for the pressure update and finally, the interlayer viscosity effect, are described in Sections 4.9 and 4.10, respectively. Section 5is devoted to the presentation of numerical examples. To test the accuracy of the scheme we consider in Example 1 (Sect. 5.1) the exact soliton solution of the LDNH0soliton (described in Appendix A). Examples 2–6 are motivated by selected experiments conducted by Iimura and Tanaka [3]. They are solved in Section 5.2 by the single-layer LDNH0test model. Examples 7 and 8 are related test cases but with limited tree height, and Examples 9–13 consider trees with properties that gradually vary with height. These cases are solved by the multilayer LDNH0model in Sections 5.3 and 5.4. Finally, in Section 6we examine possible directions for future research and address the results of the study. The practical implications of the model for environmental conservation, disaster resilience, and coastal design are discussed in this work. 1116 RAIMUND B¨ URGER ET AL. 2. Preliminaries 2.1. Basic method Let us consider a uniform mesh of cells 𝐼𝑖= [𝑥𝑖−1/2, 𝑥𝑖+1/2], where 𝑥𝑖=𝑖∆𝑥,𝑖∈Z, and time steps 𝑡𝜈=𝜈∆𝑡,𝜈∈N0. Then a first-order finite volume discretization of (1.1) can be written as 𝑊𝜈+1 𝑖=𝑊𝜈 𝑖−∆𝑡 ∆𝑥(︁𝐷𝜈,+ 𝑖−1/2+𝐷𝜈,− 𝑖+1/2)︁, 𝑖 ∈Z, 𝜈 ∈N0(2.1) (see [11] for details), where we define the numerical flux vectors 𝐷𝜈,± 𝑖+1/2:=1 2(︃𝐹(︀𝑊𝜈 𝑖+1)︀−𝐹(𝑊𝜈 𝑖) + ℬ𝜈 𝑖+1/2(︀𝑊𝜈 𝑖+1 −𝑊𝜈 𝑖)︀−(𝜎𝑖+1 −𝜎𝑖)𝒮𝜈 𝑖+1/2 ±𝑄𝜈 𝑖+1/2(︂𝑊𝜈 𝑖+1 −𝑊𝜈 𝑖−(𝜎𝑖+1 −𝜎𝑖)(︁𝒜𝜈 𝑖+1/2)︁−1𝒮𝜈 𝑖+1/2)︂)︃, (2.2) where 𝑄𝜈 𝑖+1/2is a numerical viscosity matrix, ℬ𝜈 𝑖+1/2and 𝒜𝜈 𝑖+1/2are intermediate matrices, and 𝒮𝜈 𝑖+1/2is the intermediate vector, corresponding to the states 𝑊𝜈,− 𝑖+1/2and 𝑊𝜈,+ 𝑖+1/2of 𝐵,𝐴, and 𝑆, respectively. The Jacobian matrix of the system (1.1), 𝐴:=𝜕𝐹(𝑊) 𝜕𝑊+𝐵(𝑊), encodes the linear relationship between the derivative terms and the unknowns in the system. Notice that the scheme (2.1) is not conservative. A standard choice of the viscosity matrix is the one that corresponds to the HLL flux [11] given by 𝑄(𝐴) = 𝛼0𝐼+𝛼1𝐴with 𝛼0:=𝑆R|𝑆L|−𝑆L|𝑆R| 𝑆R−𝑆L , 𝛼1:=|𝑆R|−|𝑆L| 𝑆R−𝑆L , where 𝐼is the 𝒩 ×𝒩 identity matrix. If we assume that the eigenvalues 𝜆1,𝑖+1/2, . . . , 𝜆𝒩,𝑖+1/2of 𝐴are real, then a possible choice for 𝑆Land 𝑆Ris 𝑆L= min{︀𝜆1,𝑖+1/2, . . . , 𝜆𝒩,𝑖+1/2}︀, 𝑆R= max{︀𝜆1,𝑖+1/2, . . . , 𝜆𝒩,𝑖+1/2}︀. In general 𝑆Rand 𝑆Lrepresent upper and lower bounds of the region in which the eigenvalues of the system are located. 2.2. Dry/wet fronts The term “dry/wet front” (DWF) is frequently used to address the interface between a region with fluid and a region without fluid. In the case of a tsunami simulation, the DWF describes the interface between the advancing tsunami wave and the dry land. As the wave approaches the coastline, the DWF moves inland, with the speed and behavior of the front affected by a variety of factors, such as the topography and bathymetry of the coastline, and the magnitude and duration of the wave. Accurately modeling the behavior of the DWF is evidently important for predicting the behavior of waves in coastal areas, and for assessing the risk and impact of tsunami events. A correct implementation of this front is crucial to preserve the well-balancing properties of the numerical method, this is, the method should preserve the equilibrium state of the fluid flow, where the water level is constant and the velocity is zero, in the presence of non-uniform bottom topography or other sources of external forces, which is critical for predicting the behavior of waves in coastal areas. In particular, well-balanced methods ensure that the numerical solution accurately captures the location of the DWF and the speed of its movement. Numerical methods that are not well-balanced can produce spurious oscillations or artificial numerical diffusion at the DWF, which can lead to inaccurate simulation results. Such methods can also violate the conservation laws leading to nonphysical solutions. MULTILAYER SHALLOW WATER MODEL FOR TSUNAMIS AND COASTAL FOREST INTERACTION 1117 Figure 1. Schematic (a) of the physical system, (b) of the physical system with dry fronts. 2.3. Well-balanced property Consider a simple shallow-water system of equations 𝜕𝑡ℎ+𝜕𝑥𝑞𝑢= 0, 𝜕𝑡𝑞𝑢+𝜕𝑥(︂𝑞2 𝑢 ℎ+1 2𝑔ℎ2)︂=−𝑔ℎ𝑧′ b(𝑥),(2.3) where ℎis the water depth, 𝑢is the water velocity, 𝑞𝑢:=ℎ𝑢 and 𝑧b=𝑧b(𝑥) is the channel bottom with 𝑥∈[0, 𝐿] and 𝑡∈[0, 𝑇], see Figure 1a. To derive this model, and in general any other Boussinesq system, from the Navier–Stokes equations it is assumed that ℎ > 0. However, this assumption does not allow us, for example, to simulate waves reaching a coast. One needs to extend this model to handle physically correctly situations when ℎ= 0. To this end, the definition 𝜂(𝑥, 𝑡):=ℎ(𝑥, 𝑡) + 𝑧b(𝑥) is very useful since in a steady state, 𝜂should be a constant. In fact, for 𝑞𝑢= 0 the shallow water system (2.3) reduces to 𝜕𝑡ℎ= 0, 𝑔ℎ𝜕𝑥ℎ=−𝑔ℎ𝑧′ b⇒ℎ=ℎ(𝑥), 𝜕𝑥(ℎ+𝑧b)=0. To achieve that the numerical scheme does not introduce non-physical oscillations, Castro et al. [30] introduced the following well-balanced condition called “conservation property” or “C-property”: Definition 2.1 (C-property).A numerical scheme is said to possess the C-property if it exactly reproduces the steady-state solutions 𝑞𝑢≡0, ℎ≡𝜂−𝑧b, where 𝜂is a constant such that 𝜂 > max𝑥∈[0,𝐿]𝑧b(𝑥). Consequently, for model with ℎ > 0 a stable numerical scheme should have the C-property to avoid nonphysical oscillations. However, this property does not handle dry fronts. In fact, in the situation of Figure 1b there are regions where ℎ(𝑥) is not well defined, and even if we set ℎ(𝑥) = 0 in these regions, the steady state is not satisfied because 𝜂in a dry region is higher than in a wet region. To include these cases, Castro et al. defined in [30] the following “extended C-property”: Definition 2.2 (Extended C-property).A numerical scheme is said to have the extended C-property if it exactly reproduces the steady-state solutions 𝑞𝑢≡0, ℎ(𝑥) = {︃𝜂−𝑏(𝑥) if 𝜂 > 𝑧b(𝑥), 0 otherwise. Since ℎcan be zero, we need to renormalize this quantity to avoid dividing by zero (for example in the term 𝑞2 𝑢/ℎ in the second equation of (2.3)). In general, for a given variable 𝑞divided by ℎthe corresponding quotient is approximated by 𝑞 ℎ≈√2𝑞ℎ √︀ℎ4+ max{ℎ4, 𝜀4}, 1118 RAIMUND B¨ URGER ET AL. where 0 < 𝜀 ≪1 is a small constant (in general it is chosen as 𝜀= 10−6or smaller), cf., e.g., [31]. This property is fundamental to correctly describe coastal scenarios. However, shallow water models do in general not have this property but, for the model to be compatible with the extended C-property it must be treated. This treatment will be presented below. 2.4. Hydrostatic reconstruction To ensure that our numerical scheme satisfies the extended C-property, we reformulate the shallow-water scheme to properly handle the existence of dry fronts. The main challenge is the preservation of steady states while balancing the flux and source terms, especially when the water height disappears. We proceed by rewriting system (2.3) as 𝜕𝑡𝑊+𝜕𝑥𝐹(𝑊) = 𝑆𝑧′ b, where we define the vectors 𝑊:=(︂ℎ 𝑞𝑢)︂,𝐹(𝑊):=⎛ ⎝ 𝑞𝑢 𝑞2 𝑢 ℎ+𝑔ℎ2 2⎞ ⎠,𝑆:=(︂0 𝑔ℎ)︂ and the matrix of the system, that is, the Jacobian matrix of 𝐹(𝑊), 𝐴(𝑊) = 𝜕𝐹(𝑊) 𝜕𝑊=⎡ ⎣ 0 1 𝑔ℎ −𝑞2 𝑢 ℎ22𝑞𝑢 ℎ⎤ ⎦ that has the eigenvalues 𝜆1=𝑢−√︀𝑔ℎ, 𝜆2=𝑢+√︀𝑔ℎ. The 𝐹and 𝑆-terms of the numerical flux (2.2) are given by 𝐹(𝑊𝑖+1)−𝐹(𝑊𝑖)−1 2(𝜎𝑖+1 −𝜎𝑖)(𝑆𝑖+1 +𝑆𝑖) =⎛ ⎜ ⎝ 𝑞𝜈 𝑢,𝑖+1 −𝑞𝜈 𝑢,𝑖 (︀𝑞𝜈 𝑢,𝑖+1)︀2 ℎ𝜈 𝑖 +𝑔(︀ℎ𝜈 𝑖+1)︀2 2−(︀𝑞𝜈 𝑢,𝑖)︀2 ℎ𝜈 𝑖−1−𝑔(ℎ𝜈 𝑖)2 2 ⎞ ⎟ ⎠+(︃0 𝑔 2(︀ℎ𝜈 𝑖+1 +ℎ𝜈 𝑖)︀(𝑧b,𝑖+1 −𝑧b,𝑖))︃. (2.4) At steady state, 𝑞𝜈 𝑢,𝑖 = 0 for all 𝑖, and the second component of (2.4) becomes 𝑔 2(︁(︀ℎ𝜈 𝑖+1)︀2−(ℎ𝜈 𝑖)2)︁+𝑔 2(︀(︀ℎ𝜈 𝑖+1 +ℎ𝜈 𝑖)︀(𝑧b,𝑖+1 −𝑧b,𝑖))︀ =𝑔 2(︀ℎ𝜈 𝑖+1 −ℎ𝜈 𝑖)︀(︀ℎ𝜈 𝑖+1 +ℎ𝑛 𝑖)︀+𝑔 2(︀(︀ℎ𝜈 𝑖+1 +ℎ𝜈 𝑖)︀(𝑧b,𝑖+1 −𝑧b,𝑖))︀ =𝑔 2(︀ℎ𝜈 𝑖+1 +ℎ𝜈 𝑖)︀(︀(︀ℎ𝜈 𝑖+1 +𝑧b,𝑖+1)︀−(ℎ𝜈 𝑖+𝑧b,𝑖))︀. (2.5) In absence of dry fronts this last expression is zero because at steady state, ℎ+𝑧bis constant, but in a system where ℎ= 0 is allowed, it may occur that, for example, ℎ𝜈 𝑖>0, ℎ𝜈 𝑖+1 = 0, and ℎ𝑖< 𝑧𝑖+1 −𝑧𝑖, as is shown in Figure 2. In this case the last expression in (2.5) reduces to 𝑔 2ℎ𝜈 𝑖(𝑧b,𝑖+1 −(ℎ𝜈 𝑖+𝑧b,𝑖)) <0; in other words, we will get non-physical velocities that can break our simulation. The hydrostatic reconstruction method, originally introduced by Audusse et al. [32], provides a way to handle this issue by ensuring MULTILAYER SHALLOW WATER MODEL FOR TSUNAMIS AND COASTAL FOREST INTERACTION 1119 Figure 2. Schematic of a dry front. well-balanced properties in numerical schemes ensuring that the scheme satisfies the extended C-property. The hydrostatic reconstruction is defined as ℎ− 𝑖+1/2:= max{𝑧b,𝑖 +ℎ𝑖−𝑧*,0}, ℎ+ 𝑖+1/2:= max{𝑧b,𝑖+1 +ℎ𝑖+1 −𝑧*,0},(2.6) where 𝑧*:= max{𝑧b,𝑖, 𝑧b,𝑖+1}, along with (ℎ𝑢)− 𝑖+1/2:=𝑢𝑖ℎ− 𝑖+1/2and (ℎ𝑢)+ 𝑖+1/2:=𝑢𝑖+1ℎ+ 𝑖+1/2.(2.7) Using this reconstruction we may show that ˜ 𝐺𝑖+1/2:=1 ∆𝑥(︂1 2𝑔(︁ℎ+ 𝑖+1/2)︁2−1 2𝑔(︁ℎ− 𝑖+1/2)︁2)︂ is a first-order approximation of 𝑔ℎ𝜕𝑥(ℎ+𝑧b); indeed, ˜ 𝐺𝑖+1/2=𝑔 2∆𝑥(︁ℎ+ 𝑖+1/2+ℎ− 𝑖+1/2)︁(︁ℎ+ 𝑖+1/2−ℎ− 𝑖+1/2)︁ =𝑔 2∆𝑥(ℎ𝑖+ℎ𝑖+1 +𝒪(∆𝑥))(𝜕𝑥(ℎ+𝑧b) + 𝒪(∆𝑥)), (2.8) hence this last discretization is consistent with the system and independent of the bottom function. Furthermore, by using the reconstruction (2.6), (2.7) we can see how the problem of Figure 2is now solved. For that case, 𝑧*= max{𝑧b,𝑖, 𝑧b,𝑖+1}=𝑧b,𝑖+1, hence ℎ− 𝑖+1/2= max{𝑧b,𝑖 +ℎ𝑖−𝑧b,𝑖+1,0}= 0, ℎ+ 𝑖+1/2= max{𝑧b,𝑖+1 +ℎ𝑖+1 −𝑧b,𝑖+1,0}= 0, and therefore ˜ 𝐺𝑖+1/2= 0, which demonstrates that the extended C-property is satisfied. Finally, as stated before, by using this reconstruction we may treat the system (2.3) as the (shallow water) system of conservation laws 𝜕𝑡ℎ+𝜕𝑥(ℎ𝑢)=0, 𝜕𝑡(ℎ𝑢) + 𝜕𝑥(︂ℎ𝑢2+1 2𝑔ℎ2)︂= 0. In absence of an explicit bottom term the scheme (2.1), (2.2) now simplifies to (2.1) along with 𝐷𝜈,± 𝑖+1/2=1 2(︁𝐹(︀𝑊𝜈 𝑖+1)︀−𝐹(𝑊𝜈 𝑖)±𝑄𝜈 𝑖+1/2(︀𝑊𝜈 𝑖+1 −𝑊𝜈 𝑖)︀)︁. 2.5. Linearized discontinuous non-hydrostatic model (LDNH0model) The LDNH0model of computational fluid dynamics (CFD) accounts for non-hydrostatic forces. In the present model a vertically constant profile for the horizontal velocity 𝑢and the vertical velocity 𝑤are assumed along with a linear vertical profile for the non-hydrostatic pressure 𝑝. This model is given by 𝜕𝑡ℎ+𝜕𝑥(ℎ𝑢)=0,(2.9a) 1120 RAIMUND B¨ URGER ET AL. 𝜕𝑡(ℎ𝑢) + 𝜕𝑥(︂ℎ𝑢2+1 2𝑔ℎ2+ℎ𝑝)︂=−(𝑔ℎ + 2𝑝)𝑧′ b,(2.9b) 𝜕𝑡(ℎ𝑤) + 𝜕𝑥(ℎ𝑢𝑤)=2𝑝, (2.9c) 𝜕𝑥𝑢+ 2𝑤−𝑢𝑧′ b ℎ= 0.(2.9d) Since no evolution equation exists for 𝑝, this system is solved numerically by a projection method consisting of two steps per time step of length ∆𝑡. As is stated in Sections 1.1 and 1.3, this procedure mimics Chorin’s well-known method for solving the incompressible Navier–Stokes equations. (1) In the first step we solve the system 𝜕𝑡ℎ+𝜕𝑥(ℎ𝑢)=0, 𝜕𝑡(ℎ𝑢) + 𝜕𝑥(︂ℎ𝑢2+1 2𝑔ℎ2)︂=−(𝑔ℎ)𝑧′ b, 𝜕𝑡(ℎ𝑤) + 𝜕𝑥(ℎ𝑢𝑤)=0 (2.10) that arises from (2.9) by setting the non-hydrostatic pressure 𝑝to zero and omitting the last equation (2.9d). The homogeneous version of system (2.10) can be written as a first-order system 𝜕𝑡𝑊+𝜕𝑥𝐹(𝑊) = 0,𝑊= (ℎ, ℎ𝑢, ℎ𝑤)T,𝐹(𝑊) = (︂ℎ𝑢, ℎ𝑢2+1 2𝑔ℎ2, ℎ𝑢𝑤)︂T , where the eigenvalues of the flux Jacobian matrix 𝜕𝐹(𝑊)/𝜕𝑊are given by 𝜆1=𝑢−√︀𝑔ℎ, 𝜆2=𝑢, 𝜆3=𝑢+√︀𝑔ℎ, (2.11) and which is therefore hyperbolic. The system (2.10) has the form (1.1), so we can apply the aforementioned discretization and the hydrostatic reconstruction. The intermediate updated variables are given by 𝑊𝜈+1/2 𝑖=𝑊𝜈 𝑖−∆𝑡 ∆𝑥(︁𝐷𝜈,+ 𝑖−1/2+𝐷𝜈,− 𝑖+1/2)︁. In particular, since we are using an HLL viscosity matrix, we can rewrite 𝐷𝜈,± 𝑖+1/2as 𝐷𝜈,± 𝑖+1/2=1 2(︀(1 ±𝛼1)(︀𝐹(︀𝑊𝜈 𝑖+1)︀−𝐹(𝑊𝜈 𝑖))︀±𝛼0(︀𝑊𝜈 𝑖+1 −𝑊𝜈 𝑖)︀)︀.(2.12) (2) Once the first step is calculated, we can proceed with the second step, namely the projection step. In semi-discrete form this step can be written as 1 ∆𝑡(︁𝑊𝜈+1 −𝑊𝜈+1/2)︁+(︀(∇𝑃)𝜈+1)︀T=0,(2.13) where we define ∇𝑃:= (0, 𝜕𝑥(ℎ𝑝)+2𝑝𝑧′ b,−2𝑝) = (︁0, ∇𝑃)︁, ∇𝑃:= (𝜕𝑥(ℎ𝑝)+2𝑝𝑧′ b,−2𝑝).(2.14) Next, we define 𝑋:=(︂𝑢 𝑤)︂,such that 𝑊=(︂ℎ ℎ𝑋)︂. MULTILAYER SHALLOW WATER MODEL FOR TSUNAMIS AND COASTAL FOREST INTERACTION 1121 Now (2.13) can be written as ℎ𝜈+1 −ℎ𝜈+1/2 ∆𝑡= 0,(2.15a) 1 ∆𝑡(︁(ℎ𝑋)𝜈+1 −(ℎ𝑋)𝜈+1/2)︁+(︁ ∇𝑃𝜈+1)︁T =0.(2.15b) From (2.15a) we get ℎ𝜈+1 =ℎ𝜈+1/2, whereas to solve (2.15b) in order to obtain (ℎ𝑋)𝜈+1 (that is, 𝑋𝜈+1 since ℎ𝜈+1 is known at that moment), we first need to use the constraint (2.9d), evaluated at time 𝑡𝜈+1, to obtain 𝑝𝜈+1, and thereby  ∇𝑃𝜈+1 (cf. (2.14)). Writing (2.15b) component-wise, we get ℎ𝜈+1𝑢𝜈+1 = (ℎ𝑢)𝜈+1/2−∆𝑡(︀𝜕𝑥(︀ℎ𝜈+1𝑝𝜈+1)︀+ 2𝑝𝜈+1𝑧′ b)︀,(2.16) ℎ𝜈+1𝑤𝜈+1 = (ℎ𝑤)𝜈+1/2+ 2∆𝑡𝑝𝜈+1.(2.17) On the other hand, (2.9d) can be written as 2ℎ𝑤 −ℎ𝑢(𝜕𝑥ℎ+ 2𝑧′ 𝑏) + ℎ𝜕𝑥(ℎ𝑢) = 0.(2.18) Then, replacing ℎ𝜈+1, (ℎ𝑢)𝜈+1 and (ℎ𝑤)𝜈+1 from (2.16) and (2.17) in (2.18) yields the desired equation for 𝑝𝜈+1. To state it, for sake of simplicity, we rename the time index 𝜈+ 1/2 by *. The result is 2(ℎ𝑤)*−(ℎ𝑢)*(𝜕𝑥ℎ*+ 2𝑧′ b) + ℎ*𝜕𝑥(ℎ𝑢)*+ ∆𝑡(︁𝑝𝜈+1(4 + 2𝑧′ b(𝜕𝑥ℎ*+ 2𝑧′ b)) +𝜕𝑥(︀ℎ*𝑝𝜈+1)︀(𝜕𝑥ℎ*+ 2𝑧′ b)−2ℎ*𝜕𝑥(︀𝑧′ b𝑝𝜈+1)︀−ℎ*𝜕𝑥𝑥(︀ℎ*𝑝𝜈+1)︀)︁= 0. (2.19) To solve (2.19) numerically for 𝑝𝜈+1, we discretize this equation on a uniform spatial grid 𝑥𝑖=𝑖∆𝑥, 𝑖= 1,2, . . . , 𝑁. The first-order derivatives are approximated by 𝜕𝑥(︀ℎ*𝑝𝜈+1)︀≈ℎ* 𝑖+1𝑝𝜈+1 𝑖+1 −ℎ* 𝑖−1𝑝𝜈+1 𝑖−1 2∆𝑥(2.20) and the second-order derivatives by 𝜕𝑥𝑥(︀ℎ*𝑝𝜈+1)︀≈ℎ* 𝑖+1𝑝𝜈+1 𝑖+1 −2ℎ* 𝑖𝑝𝜈+1 𝑖+ℎ* 𝑖−1𝑝𝜈+1 𝑖−1 ∆𝑥2·(2.21) Rewriting (2.19) in a discrete form in terms of unknowns 𝑝𝜈+1 𝑖, we obtain a system 𝑇𝑖,𝑖−1𝑝𝜈+1 𝑖−1+𝑇𝑖,𝑖𝑝𝜈+1 𝑖+𝑇𝑖,𝑖+1𝑝𝜈+1 𝑖+1 =𝑃0,𝑖, 𝑖 = 1, . . . , 𝑁 −1, where the coefficients 𝑇𝑖,𝑖−1,𝑇𝑖,𝑖, and 𝑇𝑖,𝑖+1 arise from the finite-difference approximations (2.20) and (2.21). In matrix form the final system can be written as 𝑇 𝑃 =𝑃0,(2.22) where 𝑃=(︀𝑝𝜈+1 0, 𝑝𝜈+1 1, . . . , 𝑝𝜈+1 𝑁)︀T(2.23) is the vector of unknowns, 𝑃0= (𝑝0,0, 𝑝0,1, . . . , 𝑝0,𝑁 )T(2.24) is the vector of right-hand sides defined by (2.19), that is 𝑝0,𝑖 = (2(ℎ𝑤)*−(ℎ𝑢)*(𝜕𝑥ℎ*+ 2𝑧′ b) + ℎ*𝜕𝑥(ℎ𝑢)*)𝑥=𝑥𝑖, 𝑖 = 0, . . . , 𝑁, (2.25) 1128 RAIMUND B¨ URGER ET AL. 𝑁 ∑︁ 𝛽=𝛼+1 𝑙𝛽𝜃𝛽𝜕𝑡ℎ𝛽+ 𝑁 ∑︁ 𝛽=𝛼+1 𝑙𝛽𝜕𝑥(ℎ𝛽𝑢𝛽)=Γ𝑁+1/2−Γ𝛼+1/2. If all layers have the same height at a given 𝑥-position, that is setting ℎ𝛽=ℎ/𝑁 for all 𝛽, we obtain (𝜕𝑡ℎ) 𝛼 ∑︁ 𝛽=1 𝜃𝛽 𝑁+ 𝛼 ∑︁ 𝛽=1 𝜕𝑥(ℎ𝛽𝑢𝛽)=Γ𝛼+1/2−Γ1/2,(4.1) (𝜕𝑡ℎ) 𝑁 ∑︁ 𝛽=𝛼+1 𝜃𝛽 𝑁+ 𝑁 ∑︁ 𝛽=𝛼+1 𝜕𝑥(ℎ𝛽𝑢𝛽)=Γ𝑁+1/2−Γ𝛼+1/2(4.2) along with the equivalence ¯ 𝜃=1 ℎ∫︁𝑧𝑁+1/2 𝑧1/2 𝜃(𝑧) d𝑧=1 ℎ 𝑁 ∑︁ 𝛽=1 ∫︁𝑧𝛽+1/2 𝑧𝛽−1/2 𝜃(𝑧) d𝑧= 𝑁 ∑︁ 𝛽=1 𝜃𝛽 𝑁· Then we can rewrite equation (4.2) as (𝜕𝑡ℎ)⎛ ⎝¯ 𝜃− 𝛼 ∑︁ 𝛽=1 𝜃𝛽 𝑁⎞ ⎠+ 𝑁 ∑︁ 𝛽=𝛼+1 𝜕𝑥(ℎ𝛽𝑢𝛽)=Γ𝑁+1/2−Γ𝛼+1/2.(4.3) To find the transfer terms we multiply (4.1) by ¯ 𝜃−∑︀𝛼 𝛽=1 𝜃𝛽/𝑁 and (4.3) by ∑︀𝛼 𝛽=1 𝜃𝛽/𝑁. Subtracting the results we obtain ⎛ ⎝¯ 𝜃− 𝛼 ∑︁ 𝛽=1 𝜃𝛽 𝑁⎞ ⎠ 𝛼 ∑︁ 𝛽=1 𝜕𝑥(ℎ𝛽𝑢𝛽)−⎛ ⎝ 𝛼 ∑︁ 𝛽=1 𝜃𝛽 𝑁⎞ ⎠ 𝑁 ∑︁ 𝛽=𝛼+1 𝜕𝑥(ℎ𝛽𝑢𝛽) =⎛ ⎝¯ 𝜃− 𝛼 ∑︁ 𝛽=1 𝜃𝛽 𝑁⎞ ⎠(︀Γ𝛼+1/2−Γ1/2)︀−⎛ ⎝ 𝛼 ∑︁ 𝛽=1 𝜃𝛽 𝑁⎞ ⎠(︀Γ𝑁+1/2−Γ𝛼+1/2)︀. (4.4) Since the boundary transfer terms are zero, i.e., Γ1/2= Γ𝑁+1/2= 0, (4.4) reduces to an identity that can be written as Γ𝛼+1/2= 𝑁 ∑︁ 𝛽=1 𝛾𝛼,𝛽𝜕𝑥(ℎ𝛽𝑢𝛽),where 𝛾𝛼,𝛽 :=⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 1− 𝛼 ∑︁ 𝑘=1 𝜃𝑘 𝑁¯ 𝜃for 𝛼≥𝛽, − 𝛼 ∑︁ 𝑘=1 𝜃𝑘 𝑁¯ 𝜃for 𝛼 < 𝛽. 4.5. Drag and inertia forces for the multilayer system A widely used definition of the drag force is 𝑓D=1 2𝜌𝐶D𝐴v|𝑢|𝑢, where 𝐴vis the effective “vertical” side face area of the object and 𝐶Dis the drag coefficient. Furthermore, to describe properties of individual trees, we use the height coordinate 𝜁that is measured from the ground surface, identified here with 𝑧b. For trees, Tanaka et al. [5] characterize 𝐶Dby 𝐶D=𝐶D,ref 𝑐tr𝑐le, MULTILAYER SHALLOW WATER MODEL FOR TSUNAMIS AND COASTAL FOREST INTERACTION 1129 where we denote by 𝑐tr𝑐le the vertical average of the product of the coefficients 𝑐tr(𝜁) and 𝑐le(𝜁) that represent the effect of the trunk and leaves of the trees, respectively: 𝑐tr𝑐le :=1 ℎ∫︁ℎ 0 𝑐tr(𝜁)𝑐le(𝜁) d𝜁 and 𝐴v=𝑑ℎ. Within a multilayer system the averages must be calculated for each layer 𝛼(𝛼= 1, . . . , 𝑁). This is done by defining 𝐶D,𝛼 :=𝐶D,ref (𝑐tr𝑐le)𝛼,where (𝑐tr𝑐le)𝛼:=1 ℎ𝛼∫︁𝑧𝛼+1/2−𝑧b 𝑧𝛼−1/2−𝑧b 𝑐tr(𝜁)𝑐le(𝜁) d𝜁, 𝛼 = 1, . . . , 𝑁. On the other hand, if the diameter of a tree at height 𝜁is 𝑑(𝜁), then we employ 𝐴v,𝛼 :=ℎ𝛼¯ 𝑑𝛼,where ¯ 𝑑𝛼:=1 ℎ𝛼∫︁𝑧𝛼+1/2−𝑧b 𝑧𝛼−1/2−𝑧b 𝑑(𝜁) d𝜁, 𝛼 = 1, . . . , 𝑁. (4.5) If we assume that the trees are symmetric along the 𝑧-axis the layer-specific transversal area 𝐴t,𝛼 of each tree is 𝐴t,𝛼 :=𝜋(¯ 𝑑𝛼)2/4 for 𝛼= 1, . . . , 𝑁. Summarizing all ingredients, we obtain that the drag force for one tree associated with layer 𝛼is given by 𝑓t,𝛼 :=1 2𝜌𝐶D,𝛼𝐴v,𝛼|𝑢𝛼|𝑢𝛼=1 2𝜌𝐶D,ref (𝑐tr𝑐le)𝛼¯ 𝑑𝛼|ℎ𝛼𝑢𝛼|𝑢𝛼, 𝛼 = 1, . . . , 𝑁. On the other hand, Tanaka et al. [5] calculate the effective density 𝑛tof the forest as 𝑛t=𝑑𝑛𝑐tr𝑐le 𝐴F , where 𝐴Fis the forest area and 𝑛denotes the number of trees divided by the forest length in flow direction. Consequently, the effective forest density 𝑛t,𝛼 and the corresponding porosity 𝜃𝛼for layer 𝛼are given by the respective expressions 𝑛t,𝛼 =¯ 𝑑𝛼𝑛(𝑐tr𝑐le)𝛼 𝐴F , 𝜃𝛼= 1 −𝑛t𝛼𝜋¯ 𝑑2 𝛼 4· Combining all ingredients we obtain the total forest drag force associated with layer 𝛼 𝑓D,𝛼 =𝑛t,𝛼 𝜃2 𝛼 𝑓t,𝛼 =1 2 𝑛t,𝛼 𝜃2 𝛼 𝜌𝐶D,ref (𝑐tr𝑐le)𝛼¯ 𝑑𝛼|ℎ𝛼𝑢𝛼|𝑢𝛼, 𝛼 = 1, . . . , 𝑁. Finally, assuming radial symmetry of the trees along the 𝑧-axis we obtain the inertia force 𝑓M,𝛼 =𝐶M 𝑛t,𝛼 𝜃𝛼 ℎ𝛼 𝜋¯ 𝑑2 𝛼 4𝜕𝑡(︂ℎ𝛼𝑢𝛼 ℎ𝛼)︂, 𝛼 = 1, . . . , 𝑁. 4.6. Gauckler–Manning friction and viscosity The friction force with respect to the ground 𝜏𝛼is present in the bottom layer only (𝛼= 1). Therefore, we define 𝜏𝛼:=𝑘1,𝛼|ℎ𝛼𝑢𝛼|ℎ𝛼𝑢𝛼=⎧ ⎨ ⎩ 𝑔𝑚2|𝑢𝛼|𝑢𝛼 𝜃𝛼ℎ1/3for 𝛼= 1, 0 for 𝛼= 2, . . . , 𝑁. 1130 RAIMUND B¨ URGER ET AL. The upper layers will be not affected directly by this force, but physically, the interaction with the ground is transferred to the upper layers by the viscosity of the fluid. This phenomenon is modeled by additional viscosity terms 𝐾𝛼+1/2and 𝐾𝛼−1/2defined by 𝐾𝛼+1/2:=−𝜂0 𝑢𝛼+1 −𝑢𝛼 ℎ𝛼+1 +ℎ𝛼 ,(4.6) where we take into account that 𝐾𝑁+1/2= 0 and 𝐾1/2=−𝜏1and 𝜂0is the viscosity constant of the fluid. These viscosity therms are a simplified expression of the general case when each layer can have different size and fluids, see [33] for further explanations. 4.7. Forces in the model Now, combining all the forces with the multilayer model, we arrive at the following governing equations of the multilayer approach: 𝜕𝑡ℎ𝛼+𝜕𝑥(ℎ𝛼𝑢𝛼) 𝜃𝛼−Γ𝛼+1/2 𝜃𝛼 +Γ𝛼−1/2 𝜃𝛼 = 0, 𝜕𝑡(ℎ𝛼𝑢𝛼) + 𝜕𝑥(︀ℎ𝛼𝑢2 𝛼)︀ 𝜃𝛼 +𝜃𝛼𝜕𝑥(ℎ𝛼𝑞𝛼)−𝜃𝛼𝜕𝑥𝑧𝛼+1/2𝑞𝛼+1/2+𝜃𝛼𝜕𝑥𝑧𝛼−1/2𝑞𝛼−1/2 =−𝑔ℎ𝛼𝜃𝛼𝜕𝑥𝜂+˜𝑢𝛼+1/2Γ𝛼+1/2 𝜃𝛼−˜𝑢𝛼−1/2Γ𝛼−1/2 𝜃𝛼−𝑓D𝛼−𝑓M𝛼+𝐾𝛼−1/2 𝜃𝛼−𝐾𝛼+1/2 𝜃𝛼 , 𝜕𝑡(ℎ𝛼𝑤𝛼) + 𝜕𝑥(ℎ𝛼𝑢𝛼𝑤𝛼) 𝜃𝛼 +𝑞𝛼+1/2−𝑞𝛼−1/2=˜𝑤𝛼+1/2Γ𝛼+1/2 𝜃𝛼−˜𝑤𝛼−1/2Γ𝛼−1/2 𝜃𝛼 , (4.7) where 𝛼= 1, . . . , 𝑁. Based on our previous calculations, we obtain the same kind of system of equations as (2.26) but this time for each layer. 4.8. Solving the first-order system Initially, our focus is on solving the first-order equations of the problem. Therefore, we assume that all pressures in (4.7) are zero during this stage of the analysis. The impact of viscosity will be incorporated in the final stages of the calculations and, hence, is disregarded in this section. Similarly to our previous approach, we derive the flow rate equation by incorporating the relevant external forces, which can be expressed as 𝜕𝑡(ℎ𝛼𝑢𝛼) + 𝜕𝑥(︀ℎ𝛼𝑢2 𝛼)︀ 𝜃𝛼−˜𝑢𝛼+1/2Γ𝛼+1/2 𝜃𝛼 +˜𝑢𝛼−1/2Γ𝛼−1/2 𝜃𝛼 =−𝑔ℎ𝛼𝜃𝛼𝜕𝑥𝜂−𝑘2,𝛼(ℎ𝛼𝑢𝛼)|ℎ𝛼𝑢𝛼|−𝑘3,𝛼(︂𝜕𝑡(ℎ𝛼𝑢𝛼)−ℎ𝛼𝑢𝛼𝜕𝑡ℎ𝛼 ℎ𝛼)︂, (1 + 𝑘3,𝛼)𝜕𝑡(ℎ𝛼𝑢𝛼)−𝑘3,𝛼𝑢𝛼𝜕𝑡ℎ𝛼+𝜕𝑥(︀ℎ𝛼𝑢2 𝛼)︀ 𝜃𝛼−˜𝑢𝛼+1/2Γ𝛼+1/2 𝜃𝛼 +˜𝑢𝛼−1/2Γ𝛼−1/2 𝜃𝛼 =−𝑔ℎ𝛼𝜃𝛼𝜕𝑥𝜂−𝑘2,𝛼(ℎ𝛼𝑢𝛼)|ℎ𝛼𝑢𝛼|, where we define the constants 𝑘2,𝛼 :=𝐶D,𝛼𝑑𝛼𝑛t,𝛼 2𝜃ℎ𝛼 and 𝑘3,𝛼 =𝐶M 𝑛t,𝛼𝜋¯ 𝑑2 𝛼 4· For every 𝛼= 1, . . . , 𝑁 we need to solve the system of balance equations ℳ𝛼𝜕𝑡⎛ ⎝ ℎ𝛼 ℎ𝛼𝑢𝛼 ℎ𝛼𝑤𝛼 ⎞ ⎠+1 𝜃𝛼 𝜕𝑥⎛ ⎜ ⎝ ℎ𝛼𝑢𝛼 ℎ𝛼𝑢2 𝛼 ℎ𝛼𝑢𝛼𝑤𝛼 ⎞ ⎟ ⎠+1 𝜃𝛼⎡ ⎣ 𝐵ℎ,𝛼 𝐵ℎ𝑢,𝛼 𝐵ℎ𝑤,𝛼 ⎤ ⎦𝜕𝑥𝑊=−⎛ ⎝ 0 𝑘2,𝛼(ℎ𝛼𝑢𝛼)|ℎ𝛼𝑢𝛼| 0⎞ ⎠,(4.8) MULTILAYER SHALLOW WATER MODEL FOR TSUNAMIS AND COASTAL FOREST INTERACTION 1131 where we define ℳ𝛼:=⎡ ⎣ 1 0 0 −𝑢𝛼𝑘3,𝛼 1 + 𝑘3,𝛼 0 0 0 1⎤ ⎦ and 𝐵ℎ,𝛼,𝐵ℎ𝑢,𝛼 and 𝐵ℎ𝑤,𝛼 represent the corresponding rows of 𝐵𝛼. Since the transfer terms depend on the velocities of all layers, 𝑊is a vector defined as 𝑊:= (ℎ1, ℎ1𝑢1, ℎ1𝑤1, ℎ2, ℎ2𝑢2, ℎ2𝑤2, . . . , ℎ𝑁, ℎ𝑁𝑢𝑁, ℎ𝑁𝑤𝑁)T∈R3𝑁. Multiplying (4.8) from the left by 𝒞𝛼:=ℳ−1 𝛼=⎡ ⎣ 1 0 0 𝑢𝛼𝑘3,𝛼/(1 + 𝑘3,𝛼) 1/(1 + 𝑘3,𝛼) 0 0 0 1⎤ ⎦ we obtain 𝜕𝑡⎛ ⎝ ℎ𝛼 ℎ𝛼𝑢𝛼 ℎ𝛼𝑤𝛼 ⎞ ⎠+1 𝜃𝛼𝒞𝛼𝜕𝑥⎛ ⎜ ⎝ ℎ𝛼𝑢𝛼 ℎ𝛼𝑢2 𝛼 ℎ𝛼𝑢𝛼𝑤𝛼 ⎞ ⎟ ⎠+1 𝜃𝛼𝒞𝛼⎡ ⎣ 𝐵ℎ,𝛼 𝐵ℎ𝑢,𝛼 𝐵ℎ𝑤,𝛼 ⎤ ⎦𝜕𝑥𝑊=−𝒞𝛼⎛ ⎝ 0 𝑘2,𝛼(ℎ𝛼𝑢𝛼)|ℎ𝛼𝑢𝛼| 0⎞ ⎠.(4.9) Now we can use the FV method described before to numerically solve the system of PDEs (4.9), but to improve the stability of the numerical scheme we discretize the remaining external forces in a semi-implicit form, this means that ℎ𝛼𝑢𝛼|ℎ𝛼𝑢𝛼|is evaluated as ℎ𝛼𝑢𝛼|ℎ𝛼𝑢𝛼| ≈ (ℎ𝛼𝑢𝛼)𝜈+1|ℎ𝛼𝑢𝛼|𝜈. Consequently, the hyperbolic scheme for layer 𝛼becomes 𝑊𝜈+1 𝛼,𝑖 =𝑁𝛼,𝑖(︂𝑊𝜈 𝛼,𝑖 −∆𝑡 ∆𝑥(︁𝐷𝜈,+ 𝛼,𝑖−1/2+𝐷𝜈,− 𝛼,𝑖+1/2)︁)︂, where the diagonal matrix 𝑁𝛼is defined analogously to (3.6); namely we here get 𝑁𝛼,𝑖 := diag(︃1,1⧸︃(︃1 + 𝑘𝜈 2,𝛼ℎ𝜈 𝛼,𝑖𝑢𝜈 𝛼,𝑖∆𝑡 1 + 𝑘3,𝛼 )︃,1)︃= diag(︃1,1 + 𝑘3,𝛼 1 + 𝑘3,𝛼 +𝑘𝜈 2,𝛼ℎ𝜈 𝛼,𝑖𝑢𝜈 𝛼,𝑖∆𝑡,1)︃, and 𝐷𝜈,± 𝛼,𝑖+1/2=1 2𝜃𝜈 𝛼,𝑖+1/2 𝒞𝜈 𝛼,𝑖+1/2(︀𝐹(︀𝑊𝜈 𝛼,𝑖+1)︀−𝐹(︀𝑊𝜈 𝛼,𝑖)︀+𝐵𝜈 𝛼,𝑖+1/2(︀𝑊𝜈 𝛼,𝑖+1 −𝑊𝜈 𝛼,𝑖)︀ ±𝑄𝜈 𝛼,𝑖+1/2(︀𝑊𝜈 𝛼,𝑖+1 −𝑊𝜈 𝛼,𝑖)︀)︀. 4.9. Projection matrix Similarly to (3.7) the elliptic problems for each layer can be written as ℎ𝜈+1 =ℎ*, (ℎ𝛼𝑢𝛼)𝜈+1 = (ℎ𝛼𝑢𝛼)*+ ∆𝑡𝑓* 𝛼(︁𝜕𝑥𝑧* 𝛼+1/2𝑞𝜈+1 𝛼+1/2−𝜕𝑥𝑧* 𝛼−1/2𝑞𝜈+1 𝛼−1/2−ℎ* 𝛼𝑞𝜈+1 𝛼)︁, (ℎ𝛼𝑤𝛼)𝜈+1 = (ℎ𝛼𝑤𝛼)*−∆𝑡(︁𝑞𝜈+1 𝛼+1/2−𝑞𝜈+1 𝛼−1/2)︁, (4.10) 1132 RAIMUND B¨ URGER ET AL. where 𝑓𝛼represents the multiplicative factors caused by the treatment of the vegetation forces, i.e., 𝑓𝛼=𝜃𝛼 1 + 𝑘3,𝛼 +𝑘2,𝛼|ℎ𝛼𝑢𝛼|∆𝑡· For 𝛼= 2, . . . , 𝑁 the constraints can be rewritten as ℎ𝛼𝑤𝛼−ℎ𝛼𝑤𝛼−1−ℎ𝛼𝑢𝛼𝜕𝑥𝑧𝛼+ℎ𝛼−1𝑢𝛼−1𝜕𝑥𝑧𝛼−1+ℎ𝛼 2𝜕𝑥(ℎ𝛼𝑢𝛼+ℎ𝛼−1𝑢𝛼−1) = 0 (4.11) and for the first layer (𝛼= 1) as ℎ1𝑤1−ℎ1𝑢1𝜕𝑥𝑧1+ℎ1 2𝜕𝑥(ℎ1𝑢1) = 0.(4.12) Typically, the variables are evaluated in volumes while the pressures are defined on the edges. To apply the projection method, we discretize the terms (ℎ𝛼𝑢𝛼)𝜈+1 and (ℎ𝛼𝑤𝛼)𝜈+1 on the edges and substitute them into the constraint equations (4.11) and (4.12). Alternatively, we can substitute these terms into (4.11) and (4.12) first and then discretize them “at 𝑖+ 1/2”, that is, on the edges. In this context, we have opted for the latter approach, which yields the equation 4𝑁2𝑃* ∆𝑡+𝑞𝛼−1/2(︀8𝑁2+𝑓𝛼𝜑2 𝛼+𝑓𝛼−1𝜑2 𝛼−1−ℎ𝜕𝑥(𝑓𝛼𝜑𝛼) + ℎ𝜕𝑥(𝑓𝛼−1𝜑𝛼−1))︀ −(︀𝜕𝑥𝑞𝛼−1/2)︀(ℎ𝜕𝑥(𝑓𝛼ℎ) + ℎ𝜕𝑥(𝑓𝛼−1ℎ)) −𝜕𝑥𝑥𝑞𝛼−1/2(︀ℎ2(𝑓𝛼+𝑓𝛼−1))︀ +𝑞𝛼−3/2(︀2𝑁2+𝑓𝛼−1𝜑2 𝛼−1+ℎ𝜕𝑥(𝑓𝛼−1𝜑𝛼−1))︀ +(︀𝜕𝑥𝑞𝛼−3/2)︀(︀ℎ𝑓𝛼−1𝜑𝛼−1+ℎ2𝜕𝑥(𝑓𝛼−1) + ℎ𝑓𝛼−1𝜑𝛼−1/2)︀+(︀𝜕𝑥𝑥𝑞𝛼−3/2)︀(︀ℎ2𝑓𝛼−1)︀ +𝑞𝛼+1/2(︀2𝑁2+𝑓𝛼𝜑2 𝛼−ℎ𝜕𝑥(𝑓𝛼𝜑𝛼))︀ +(︀𝜕𝑥𝑞𝛼+1/2)︀(︀ℎ𝑓𝛼𝜑𝛼+ℎ2𝜕𝑥(𝑓𝛼) + ℎ𝑓𝛼𝜑𝛼−1/2)︀+𝜕𝑥𝑥𝑞𝛼+1/2(︀ℎ2𝑓𝛼)︀= 0, (4.13) where 𝑃*are the independent terms. The system (4.13) is solved numerically by an iterative Jacobi method. In this way we only need to write out the part of the matrix associated with 𝑞𝛼−1/2. If we denote this matrix by 𝑀= (𝑀𝑖,𝑗) then 𝑀𝑖,𝑖 = 8𝑁2+𝑓𝛼𝜑2 𝛼+𝑓𝛼−1𝜑2 𝛼−1−ℎ𝑖+1/2𝜕𝑥(︀𝑓𝛼,𝑖+1/2𝜑𝛼,𝑖+1/2)︀+ℎ𝑖+1/2𝜕𝑥(︀𝑓𝛼−1,𝑖+1/2𝜑𝛼−1)︀ +2 ∆𝑥2ℎ2 𝑖+1/2(︀𝑓𝛼,𝑖+1/2+𝑓𝛼−1,𝑖+1/2)︀, 𝑀𝑖,𝑖+1 =−1 2∆𝑥(︀ℎ𝑖+1/2𝜕𝑥(︀𝑓𝛼,𝑖+1/2ℎ𝑖+1/2)︀+ℎ𝑖+1/2𝜕𝑥(︀𝑓𝛼−1,𝑖+1/2ℎ𝑖+1/2)︀)︀ −1 ∆𝑥2ℎ2 𝑖+1/2(︀𝑓𝛼,𝑖+1/2+𝑓𝛼−1,𝑖+1/2)︀, 𝑀𝑖,𝑖−1=1 2∆𝑥(︀ℎ𝑖+1/2𝜕𝑥(︀𝑓𝛼,𝑖+1/2ℎ𝑖+1/2)︀+ℎ𝑖+1/2𝜕𝑥(︀𝑓𝛼−1,𝑖+1/2ℎ𝑖+1/2)︀)︀ −1 ∆𝑥2ℎ2 𝑖+1/2(︀𝑓𝛼,𝑖+1/2+𝑓𝛼−1,𝑖+1/2)︀ and 𝑀𝑖,𝑗 = 0 for |𝑖−𝑗|>1. Notice that 𝜕𝑥ℎ=𝜑𝛼−𝜑𝛼−1, where 𝜑𝛼= 2𝑁𝜕𝑥(︂𝑧𝛼−1/2+ℎ𝛼 2)︂· As stated before, once the values of 𝑞𝛼−3/2,𝑞𝛼−1/2and 𝑞𝛼+1/2have been found for each layer, we must update the values of 𝑢𝛼and 𝑤𝛼for each layer using (4.10). MULTILAYER SHALLOW WATER MODEL FOR TSUNAMIS AND COASTAL FOREST INTERACTION 1133 4.10. Adding viscosity and friction At this point we have variables evaluated at time step 𝑡𝜈+1 but we must add the viscosity effect, so we will rename this time as *and because we are only correcting the velocity, thereby omitting the superscript in ℎ𝛼. Thus, we obtain ℎ𝛼𝑢𝜈+1 𝛼=ℎ𝛼𝑢* 𝛼+∆𝑡 𝜃𝛼 𝐾𝜈+1 𝛼−1/2−∆𝑡 𝜃𝛼 𝐾𝜈+1 𝛼+1/2. By using the definition (4.6) we get ℎ1𝑢𝜈+1 1=ℎ1𝑢* 1−∆𝑡𝑘1|ℎ1𝑢* 1|ℎ1𝑢𝜈+1 1+ ∆𝑡𝜂0 2𝜃𝛼 𝑢𝜈+1 2−𝑢𝜈+1 1 ℎ1 , ℎ𝛼𝑢𝜈+1 𝛼=ℎ𝛼𝑢* 𝛼−∆𝑡𝜂0 2𝜃𝛼 𝑢𝜈+1 𝛼−𝑢𝜈+1 𝛼−1 ℎ𝛼 + ∆𝑡𝜂0 2𝜃𝛼 𝑢𝜈+1 𝛼+1 −𝑢𝜈+1 𝛼 ℎ𝛼 , 𝛼 = 2, . . . , 𝑁 −1, ℎ𝑁𝑢𝜈+1 𝑁=ℎ𝑁𝑢* 𝑁−∆𝑡𝜂0 2𝜃𝛼 𝑢𝜈+1 𝑁−𝑢𝜈+1 𝑁−1 ℎ𝑁· For the unknowns 𝑢𝜈+1 1, . . . , 𝑢𝜈+1 𝑁we obtain the linear system of equations (︂1 + 𝜂0∆𝑡 2𝜃𝛼ℎ2 1 + ∆𝑡𝑘1|ℎ1𝑢* 1|)︂ℎ1𝑢𝜈+1 1−𝜂0∆𝑡 2𝜃𝛼ℎ2 𝛼 ℎ𝛼𝑢𝜈+1 2=ℎ1𝑢* 1, (︂1 + 𝜂0∆𝑡 𝜃𝛼ℎ2 𝛼)︂ℎ𝛼𝑢𝜈+1 𝛼−𝜂0∆𝑡 2𝜃𝛼ℎ2 𝛼 ℎ𝛼𝑢𝜈+1 𝛼−1−𝜂0∆𝑡 2𝜃𝛼ℎ2 𝛼 ℎ𝛼𝑢𝜈+1 𝛼+1 =ℎ𝛼𝑢* 𝛼, 𝛼 = 2, . . . , 𝑁 −1, (︂1 + 𝜂0∆𝑡 2𝜃𝛼ℎ2 𝑁)︂ℎ𝑁𝑢𝜈+1 𝑁−𝜂0∆𝑡 2𝜃𝛼ℎ2 𝛼 ℎ𝛼𝑢𝜈+1 𝑁−1=ℎ𝑁𝑢* 𝑁. This system has a tridiagonal matrix and is solved by a Thomas algorithm. 5. Numerical results 5.1. Example 1: convergence test As a first test of the accuracy of the scheme we consider the exact soliton solution of the LDNH0model, which is described in Appendix A. The computational domain is the 𝑥-interval 𝑋= [−25,25], which is subdivided into 𝒥subintervals of length ∆𝑥= 50/𝒥, and we let 𝜙(𝑥𝑖, 𝑡) denote the numerical approximation of the exact value 𝜙exact(𝑥𝑖, 𝑡), where 𝑥𝑖is the 𝑖−th discrete spatial point in the grid. We measure the 𝐿1 error in 𝜙at simulated time 𝑡= 10 s as follows: 𝑒𝜙(𝑡):=1 𝒥 𝒥 ∑︁ 𝑖=1𝜙(𝑥𝑖, 𝑡)−𝜙exact(𝑥𝑖, 𝑡); this is done for 𝜙=ℎ,𝜙=ℎ𝑢, and 𝜙=ℎ𝑤. We utilize a CFL number of 0.8 and obtain the errors displayed in Table 1. We observe that the error decreases at a rate slightly smaller than one, in agreement with the formal first-order accuracy of the numerical scheme. 5.2. Examples 2–6: Iimura–Tanaka experiments and tests for the LDNH0single-layer model The model and numerical method are motivated by a series of experiments reported by Iimura and Tanaka [3] that represents a scale 1:100 scenario for a real-world tsunami. The experimental setup consists in a channel of width 0.4 m and length 15 m with a bottom topography representing a beach (coastal area) with two slopes, see Figure 1 of [3] and our Figure 4. At the seaward end of the channel, at 𝑥= 0 m, a 1134 RAIMUND B¨ URGER ET AL. Table 1. Example 1: convergence test (comparison with an exact soliton solution). Number 𝐿1error 𝐿1error 𝐿1error of cells 𝒥𝑒ℎRate 𝑒ℎ𝑢 Rate 𝑒ℎ𝑤 Rate 50 1.17E−02 – 3.99E−02 – 1.14E−02 – 100 5.30E−03 1.14 1.78E−02 1.16 6.20E−03 0.88 200 3.00E−03 0.82 1.01E−02 0.82 3.80E−03 0.71 400 1.70E−03 0.82 5.70E−03 0.83 2.20E−03 0.79 800 9.11E−04 0.90 3.00E−03 0.93 1.20E−03 0.87 1600 4.75E−04 0.94 1.60E−03 0.91 6.56E−04 0.87 Figure 4. Schematic of the experiment, showing the positions of the nine measurement points 𝐺1, . . . , 𝐺9. The positions of 𝐺1to 𝐺6are fixed, 𝐺7and 𝐺8are in front of and behind the vegetation, and 𝐺8is in the middle of the vegetation for uniform arrangements (in this study, Cases 1, 3 and 5 of [3]) and at the boundary of two different tree densities in the combined arrangements (in this work, Cases 8 and 13 of [3]). The specific situation in this plot with 𝐺7= 10.36 m, 𝐺8= 10.36 m and 𝐺9= 11.36 m corresponds to Case 1. Notice that the vertical scale is five times larger than the horizontal. wave-making plate is located, and between 𝑥= 10.36 m and 𝑥= 11.36 m various arrangements of vertical cylinders, each with a diameter of 𝑑= 0.005 m can be placed to model the coastal vegetation (see Fig. 5). The level of water at rest is 0.4 m. The run-up height, water level, and force acting on a cylinder were measured at nine different points (𝐺1, . . . , 𝐺9; see Figure 4for the corresponding 𝑥-positions). Among the 15 different distributions of coastal vegetation tested in [3], Cases 1–15, we used five, namely Cases 1, 3, 5, 8, and 13, for numerical simulation. To specify the corresponding parameters, we recall from [3] that for given a tree distribution such as the one drawn in Figure 5, the thickness of vegetation is calculated as 𝑑𝑛=2 √3𝐷2 f 𝑊f𝑑×105+2 √3𝐷2 b 𝑊b𝑑×105,(5.1) where the factor 105adjusts a unit of 𝑑𝑛because 𝐷and 𝑊are measured in millimeters and 𝑑in meters. In the experiments, 𝑑𝑛was set to a constant 231 in all experiments [3]. The width of the channel is 0.4 m; then, the tree density, measured in trees per square metre, can be calculated as 𝑛t=𝑑𝑛 𝑊×0.4 m· The wave-making plate at 𝑥= 0 m generates a solitary wave with a height of 3.14 cm in the vicinity of the left boundary. MULTILAYER SHALLOW WATER MODEL FOR TSUNAMIS AND COASTAL FOREST INTERACTION 1135 Figure 5. Schematic of the distribution of trees in the experiment, as seen from above. The (model) forest in flow direction has the total length 𝑊, and is subdivided into a front part and a back part of lengths 𝑊fand 𝑊b, respectively. These parts may be equipped with different models of vegetation. Table 2. Parameters and drag coefficients for Cases 1, 3, 5, 8, and 13 from Iimura and Tanaka [3]. The table includes tree spacing, forest length, drag coefficient values, and tree densities for each case. Examples 2–6 correspond to the simulations in Figure 6. Example 2 Example 3 Example 4 Example 5 Example 6 Parameter Case 1 Case 3 Case 5 Case 8 Case 13 𝐷b[mm] 50 30 10 40 20 𝐷f[mm] 0 0 0 20 40 𝑊b[mm] 0 360 40 320 80 𝑊f[mm] 1000 0 0 80 320 𝑥F[m] 10.32 11.00 11.32 10.96 10.96 𝐶D(𝑥7) 0.71 0.66 1.73 0.83 0.76 𝐶D(𝑥8,f) 0.94 0.79 1.32 1.07 0.76 𝐶D(𝑥8,b) 0.94 0.79 1.32 0.84 0.92 𝐶D(𝑥9) 0.77 0.94 2.23 0.74 1.26 𝑛t[m−2] 577.5 1604.16 14437.5 𝑛t,f= 3609.4 𝑛t,b= 902.3 𝑛t,f= 902.3 𝑛t,b= 3609.4 The parameters defining the five cases considered herein are summarized in Table 2. The drag coefficient 𝐶Dis defined as a function of 𝑥obtained by piecewise linear interpolation of values of 𝐶Dspecified at 𝑥7, 𝑥8, and 𝑥9, or at 𝑥8,fand 𝑥8,binstead of 𝑥8in the case that the “front” and “back” limiting values of 𝐶D differ (see Fig. 5). On the other hand, in [34], the authors analyse and compare different approaches for describing and generating solitary waves by plates in a flume. We found that the profile that best matched the initial condition was the solution found by Rayleigh [35] that is given by 𝜂(𝑥, 0) = 𝜂0sech2(𝛽𝑥) with 𝛽=√︃3𝜂0 4ℎ2 0(ℎ0+𝜂0),(5.2) where 𝜂0is the wave amplitude and ℎ0height of water. This means that the initial wave velocity is 𝑢(𝑥, 0) = 𝑐𝜂(𝑥, 0) ℎ(𝑥, 0) with 𝑐=√︀𝑔(ℎ0+𝜂0). 1136 RAIMUND B¨ URGER ET AL. Figure 6. Examples 2–6: experimental data (red dots) and LDNH0simulations for (a) Example 2 ([3], Case 1), (b) Example 3 ([3], Case 3), (c) Example 4 ([3], Case 5), (d) Example 5 ([3], Case 8) and (e) Example 6 ([3], Case 13). In order to fit the wave into the channel we slightly extend the computational domain to 𝑥 < 0 (Fig. 7). In addition to the experimental data we include in Figure 6numerical simulations of the maximal water level. All numerical simulations for Examples 2–16 have been obtained with ∆𝑥= 0.02 m and are based on a CFL number of 0.8. The simulations of Examples 2–6 in Figure 6have been obtained by the LDNH0. For Case 3 of [3] (our Example 3) we also display in Figure 8b the simulated water height at 𝑥=𝑥9= 11.36 m, the position of 𝐺9, as a function of time. 5.3. Examples 7 and 8: Iimura–Tanaka experiment and LDNH0multi-layer tests. Having calibrated the model using the experimental data, we are now ready to simulate diverse scenarios that account for trees with properties that vary along the vertical axis. We use the same initial condition as for Examples 2–6. All simulations employ the same initial condition for the water waves, namely the one given by (5.2). Now the vegetation parameters depend on height 𝑧. Equation (5.1) is replaced by 𝑑𝑛,𝛼 =2 √3𝐷2 f 𝑊f¯ 𝑑𝛼×105+2 √3𝐷2 b 𝑊b¯ 𝑑𝛼×105, 𝛼 = 1, . . . , 𝑁, that is we assume that 𝑑𝑛is specific for each layer 𝛼. The average tree diameter ¯ 𝑑𝛼is defined in (4.5). On the other hand, 𝐷fand 𝐷bare calculated from the centers of the trees and we assume trees have layer-wise cylindrical symmetry so they do not change with respect to 𝜁. Likewise, 𝑊fand 𝑊bdo not depend on 𝜁or MULTILAYER SHALLOW WATER MODEL FOR TSUNAMIS AND COASTAL FOREST INTERACTION 1137 Figure 7. Numerical simulations of (a) Case 1, (b) Case 3, (c) Case 5 and (d) Case 13. Figure 8. Examples 2–6: (a) initial condition of the simulation, (b) Example 3 ([3], Case 3): temporal evolution at 𝑥= 11.36 m compared with experimental data. 𝑧either. Furthermore, to include the effects of branches and leaves given by 𝑐tr(𝜁)𝑐le(𝜁) (see Sect. 4.5) one usually defines 𝑑𝑛,all :=𝑑𝑛×𝑐tr𝑐le. For the multilayer model we generalize this as 𝑑𝑛,all,𝛼 =𝑑𝑛,𝛼 ×(𝑐tr𝑐le)𝛼, 𝛼 = 1, . . . , 𝑁. Next, we investigate the impact of varying tree heights. In Example 7 we assume that the trees have a height of 0.03 m (measured from ground). This information is incorporated by setting 𝑑(𝜁) = {︃0.03 m for 𝜁 < 0.03 m, 0 m for 𝜁≥0.03 m. In this case differences in the numerical solution in dependence of the number of layers do become visible in the forest areas (Example 7, see Fig. 9a). 1144 RAIMUND B¨ URGER ET AL. Inserting all these expressions into (2.9c) we end up with the equation −2𝐴 ℎ(𝜉)−2ℎ3 0𝛾2𝛽2 𝜂0ℎ2(𝜉)−2ℎ2 0𝛾2𝛽2 ℎ2(𝜉)+3ℎ2 0𝛾2𝛽2 𝜂0ℎ(𝜉)+2ℎ0𝛾2𝛽2 ℎ(𝜉)−𝛾2𝛽2ℎ(𝜉) 𝜂0 +2𝛾2 ℎ2(𝜉)+𝑔ℎ(𝜉) = 0 or equivalently, (︂𝑔−𝛾2𝛽2 𝜂0)︂ℎ3(𝜉) + (︂2ℎ0𝛾2𝛽2+3ℎ2 0𝛾2𝛽2 𝜂0 −2𝐴)︂ℎ(𝜉) + (︂2𝛾2−2ℎ3 0𝛾2𝛽2 𝜂0 −2ℎ2 0𝛾2𝛽2)︂= 0. From the requirement that the coefficients of ℎ3(𝜉), ℎ(𝜉), and the independent term should vanish independently we find that 𝛽=√︂𝜂0 ℎ2 0(ℎ0+𝜂0)and 𝑐=√︀𝑔(ℎ0+𝜂0). Writing again (𝑥, 𝑡) instead of 𝜉we obtain the soliton solution ℎ(𝑥, 𝑡) = ℎ0+𝜂0sech2(𝛽(𝑥−𝑐𝑡)), 𝑢(𝑥, 𝑡) = 𝑐(︂1−ℎ0 ℎ)︂, 𝑤(𝑥, 𝑡) = 𝑐𝛽ℎ0tanh(𝛽(𝑥−𝑐𝑡))(ℎ−ℎ0) ℎ, 𝑝(𝑥, 𝑡) = 𝑔ℎ0(3ℎ0+ 2𝜂0) 2ℎ−(ℎ0𝑐)2 ℎ2−𝑔ℎ 2, (A.1) see Figure A.1. In particular, here we showed a simple way to find this soliton because we knew the functional form of ℎ(𝜉) and it is easy to verify that this is solution of (2.9), but in order to find ℎ(𝜉) we used the sech method, and more complex models will require full use of it. Figure A.1. Graphs of LDNH0soliton solution with ℎ0= 1 and 𝜂0= 0.2.