## 1 Introduction

Coastal areas have a dynamic morphological nature driven by non-trivial interactions between sediment transport, bed morphodynamics, and water waves forced by astronomical tides, winds, and long-wave currents. Changes in coastal morphology caused by natural and anthropogenic forces have the potential to negatively affect coastal infrastructure and environment. For example, structural integrity of piers, levees and other coastal infrastructure can be compromised by excessive erosion of bed due to scouring. Moreover, sediment transport and bed morphodynamic processes play an important role in harbor planning and construction since excessive sediment deposition in a harbor may significantly increase its operating costs due to necessary dredging. Environmental concerns include shoreline and beach erosion that may damage natural habitats of endangered protected species, and effects of sediment transport on contaminants, i.e. sediment deposits may serve as dangerous contaminant sinks or sources depending on surrounding physico-chemical conditions. This evidence suggests that mathematical modeling of hydro-morphodynamic processes in coastal areas required to forecast sediment transport and bed evolution has clear engineering relevance.

Sediment transport and bed morphodynamic processes are driven by water flow properties, such as the flow velocity and turbulence, which, in turn, are affected by changes in bed surface levels. Therefore, any mathematical modeling of hydro-morphodynamic processes involves a coupling between a hydrodynamic model, that describes water waves and motion, and a sediment transport and bed evolution model, that resolves changes in bed topology driven by sediment erosion, transport, and deposition rates. A widely used variation of mathematical models for hydro-morphodynamic processes is formed by coupling the nonlinear shallow water equations with the sediment continuity Exner equation. Numerical approaches for solving the resulting set of equations can be fully coupled or decoupled, and use structured or unstructured computational grids. The use of unstructured grids can be highly advantageous as they provide the ability of local grid refinement near important bathymetric features and structures. This ability can improve accuracy while maintaining lower computational costs as compared to models that use structured grid methods. Moreover, unstructured grids have better capacity to capture irregular geometries of coastal areas, which is a great advantage over structured grids when hydro-morphodynamic processes are modeled around coastal areas.

Examples of decoupled approaches over unstructured grids with discontinuous Galerkin methods include Kubatko *et al.* kubatko_etal_2006 and Izem *et al.* izem_etal_2012. A decoupled approach suggests that the nonlinear shallow water equations and the Exner equation are solved separately from each other. In cases where the morphodynamic model has time scales much longer than the hydrodynamic model, updates in bed elevation may be done every time steps of the hydrodynamic model kubatko_etal_2006; izem_etal_2012

. Although this may provide the opportunity to reduce the amount of computational resources required to run the decoupled model, the method may not be suitable for rapidly evolving beds. In this case, a fully coupled model that solves the hydrodynamic and morphodynamic models simultaneously is more fitting. The resulting coupled model forms a system of hyperbolic nonconservative partial differential equations due to the presence of a nonconservative product in the source term. This fact adds a degree of complexity to the coupled model’s numerical solution algorithm. Among examples of discontinuous Galerkin formulations for the coupled nonconservative system are Tassi

*et al.*tassi_etal_2008, Rhebergen

*et al.*rhebergen_2008, and Mirabito

*et al.*mirabito_etal_2011. A major detail of these methods is the special treatment of the nonconservative product term developed using the theory of Dal Maso

*et al.*dal_maso_etal_1995.

The choice of the nonlinear shallow water equations is popular for a number of reasons: multiple numerical solution algorithms have been developed for these equations (e.g. discontinuous Galerkin implementations in Aizinger and Dawson aizinger_and_dawson_2002, Kubatko *et al.* kubatko_etal_2006_nswe), a track record of successful application in real world scenarios (e.g. storm surge modeling in Dawson *et al.* (dawson_etal_2011)), the ability of these equations to handle wetting-drying phenomenon that is important for coastal applications (e.g. Bunya *et al.* (bunya_etal_2009)), efficient parallelization strategies (e.g. hybrid MPI+OpenMP, and HPX parallelization in Bremer *et al.* (bremer_etal_2019)), and the ability to capture wave breaking in surf zones. Although the nonlinear shallow water equations provide this multitude of advantages, their lack of ability to capture dispersive wave effects can be a major disadvantage when water wave dynamics must be modeled in areas where wave dispersion is prevalent. An alternative depth-averaged hydrodynamic model that can capture these effects is formed by the Green-Naghdi equations developed in (green_naghdi_1976).

The capacity to capture dispersive wave effects comes, however, with a greater analytical and numerical complexity. Among numerical solution algorithms proposed for the Green-Naghdi equations, a few have been based on the Strang operator splitting technique (e.g. Bonneton *et al.* (bonneton_etal_2011), Samii and Dawson (samii_and_dawson_2018)). In this approach the Green-Naghdi equations are split into two parts: (1) the nonlinear shallow water equations, and (2) the dispersive correction part of the equations. A numerical solution operator for the Green-Naghdi equations is then defined as a successive application of numerical solution operators for these two parts. Although numerical solution algorithms for the two parts do not have to employ the same discretization method (e.g. Lannes and Marche (lannes_and_marche_2015) use a finite volume method for the first part and a finite difference method for the second part), Duran and Marche (duran_and_marche_2017) use a discontinuous Galerkin method for both parts, and Samii and Dawson (samii_and_dawson_2018) use a hybridized discontinuous Galerkin method to discretize both parts. The operator splitting approach provides a possibility to switch between the nonlinear shallow water equations and the Green-Naghdi equations when modeling water flow dynamics. The switching to the nonlinear shallow water equations can be simply done by not applying the dispersive correction part in areas where the Green-Naghdi equations provide a less accurate model, e.g. in surf zones where wave breaking occurs (bonneton_etal_2011).

This work aims to introduce a dispersive wave hydro-morphodynamic model by coupling the Green-Naghdi equations with the sediment continuity Exner equation, and to develop numerical solution algorithms for the model. Major motivation for the derivation of this model is its application in a future work to forecast morphodynamic evolution of coastal areas due to dispersive water waves. A significant portion of this work comprises the development of a massively parallel solver that uses the presented numerical solution algorithms. The solver extends a C++ software package ^{1}^{1}1The software is under development on the date of the publication, and can be accessed at www.github.com/UT-CHG/dgswemv2. Should there be any questions, comments, or suggestions, please contact the developers through the repository issues page. developed by Bremer and Kazhyken, and has the capacity to execute numerical simulations of water waves using discontinuous Galerkin discretizations of the nonlinear shallow water and Green-Naghdi equations.

This paper is organized as follows. In Section 2, the governing equations are presented for the developed mathematical model. The numerical solution algorithms, using discontinuous Galerkin methods over unstructured grids, are introduced in Section 3 both for the decoupled and coupled models. Section 4 presents a numerical test run which is used to validate the developed algorithms. The numerical setting aims to model flow and sediment transport induced by a solitary wave over a sloping beach and compare the results against the experiment conducted by Sumer *et al.* (sumer_etal_2011). Final conclusions are presented in Section 5.

## 2 Governing equations

For purposes of this work, a body of water is represented by a domain filled with water as an incompressible, homogeneous, inviscid fluid. In this description, stands for the horizontal spatial dimension that can take values 1 or 2, represents the time variable, and are the top and bottom boundaries of the domain, respectively, is the characteristic length, and is the reference depth (cf. Fig.1). It is assumed that and can be represented as graphs, and fluid particles do not cross the boundaries. Both boundaries vary with time: due to sediment transport and bed morphodynamic processes, as the evolving free surface of the body of water. The bathymetry, , and the free surface elevation, , of the body of water are used in the parameterization of and :

(1a) | ||||

(1b) |

and the domain is defined as a set of points where .

Motion of water over an erodible bed and subsequent sediment transport and bed surface evolution are highly interactive processes. Water flow parameters, such as the flow velocity and turbulence, determine the rates of sediment erosion, transport, and deposition that drive changes in bed relief; and these changes, in turn, affect the flow parameters. Therefore, any mathematical modeling of these interrelated hydro-morphodynamic processes involves some sort of coupling between a hydrodynamic model, which governs the changes in flow parameters, and a sediment transport and bed morphodynamic model, which determines the sediment erosion, transport, and deposition rates, and the subsequent changes in bed levels.

### 2.1 Hydrodynamic model

Defining the shallowness parameter , the shallow water flow regime is in action when . Under the assumption of the shallow water flow regime, the Green-Naghdi equations, a depth-averaged hydrodynamic model, provide a sufficiently accurate approximation to water flow dynamics within the domain while maintaining the ability to capture wave dispersion effects (bonneton_etal_2011). A single parameter variation of the Green-Naghdi equations introduced by Bonneton *et al.* in (bonneton_etal_2011) are defined over a horizontal domain as

(2) |

where

(3) |

is the water velocity represented by a

dimensional vector,

is the water depth represented by the mapping and assumed to be bounded from below by a positive value, comprises additional source terms for the momentum continuity equation, e.g. the Coriolis, bottom friction, and surface wind stress forces, is the acceleration due to gravity,is the identity matrix, and where the wave dispersion effects are introduced into the model through the dispersive term

(4) |

In this description, is defined through

(5) |

where operators and are

(6a) | ||||

(6b) |

with operators and defined as

(7a) | ||||

(7b) |

In this description, is a parameter that is used to optimize dispersive properties of the presented hydrodynamic model. By adjusting , the difference between the phase and group velocities coming from the Stokes linear theory and the Green-Naghdi equations can be minimized. A common strategy aims at minimizing the averaged variation over some range of wave number values (bonneton_etal_2011).

### 2.2 Sediment transport and bed morphodynamic model

Among modes of sediment transport are bed-load, suspended-load, and wash-load transport. In the presented work, the developed model is limited to bed-load transport, where sediment particles slide, roll, and saltate due to shearing forces from the surrounding fluid while staying sufficiently close to bed. The sediment continuity Exner equation provides a mathematical model that describes morphological evolution of bed due to sediment transport phenomena (exner_1925). In a morphodynamic model limited to bed-load transport, the equation states that change of in time is equal to the divergence of the sediment flux :

(8) |

where is an empirically defined function (diaz_etal_2008). Intuitively, sediment transport occurs in the flow direction; therefore,

(9) |

where is the unit flow velocity vector, and is the magnitude of the sediment flux represented by an empirical formula. A number of empirical models have been proposed for ; most of them may be represented as (see (diaz_etal_2008; cordier_etal_2011) and all references therein)

(10) |

where and is an empirical equation, e.g. the Grass model takes as a constant calibrated for the application under investigation and sets (grass_1981). There are a number of other empirical expressions for , e.g. Meyer-Peter and Mueller (meyer_and_muller_1948), Fernandez Luque and Van Beek (luque_beek_1976), Nielsen (nielsen_1992), Ribberink (ribberink_1998). The choice of the empirical representation of is judicious and influenced by the application.

## 3 Numerical methods

Discontinuous Galerkin finite element methods are used for discretizing the governing equations. This choice facilitates the use of unstructured meshes that are well suited for irregular geometries of coastal areas.

### 3.1 Notation and functional setting

The problem domain is partitioned into a finite element mesh that provides an approximation to the domain:

(11) |

where the subscript stands for the mesh parameter represented by the diameter of the smallest element in the mesh. The set of all faces of elements of the mesh, , and the set of all edges of the mesh skeleton, , are defined as

(12a) | ||||

(12b) |

Note that in the common element faces appear only once but in they are counted twice.

To develop variational formulations of the governing equations, inner products are defined for finite dimensional vectors and through:

(13a) | ||||

(13b) |

for and .

An approximating space of trial and test functions is chosen as the set of square integrable functions over such that their restriction to an element of the mesh belongs to , a space of polynomials of degree at most with support in :

(14) |

and, similarly, an approximation space over the mesh skeleton is chosen as

(15) |

### 3.2 Decoupled model

In the decoupled model method the Green-Naghdi and Exner equations are solved separately. After the flow parameters are evolved in time according to the hydrodynamic model for a number of time steps, the bed surface elevation is updated with the use of the morphodynamic model and fed back into the hydrodynamic model to continue the evolution of the flow parameters until the next bed surface elevation update. In cases where the time scales in the hydrodynamic model are much shorter than the time scales in the morphodynamic model, the bed surface elevation does not need to be updated every time step of the hydrodynamic model. In some cases the bed update may happen every time steps of the hydrodynamic model (kubatko_etal_2006; izem_etal_2012). The ability to save computational resources is the main advantage of the decoupled model method. However, this method may be unsuitable if the time scales in the hydrodynamic and morphodynamic models are comparable, e.g. in the case of a dam break (kubatko_etal_2006; izem_etal_2012).

The Green-Naghdi equations presented in Eq.(2) can be treated numerically with the use of the well-known Strang operator splitting technique (bonneton_etal_2011; samii_and_dawson_2018). The equation is split into: (1) the nonlinear shallow water equations by dropping the dispersive term of the equation, and (2) the dispersive correction part where the wave dispersion effects on flow velocities are introduced into the model through the dispersive term. If is a numerical solution operator for the nonlinear shallow water equations, i.e. propagates numerical solution by a time step , and, similarly, is a numerical solution operator for the dispersive correction part, then the second-order Strang operator splitting technique (strang_1968) states that a numerical solution operator for the Green-Naghdi equations can be approximated as

(16) |

where is a second-order temporal discretization if both and use a second-order time discretization method.

A numerical solution operator for the nonlinear shallow water equations is developed using a discontinuous Galerkin finite element formulation. Therefore, an approximate solution must satisfy the variational formulation

(17) |

where and , is a single valued approximation to over element faces, called the numerical flux, and is the unit outward normal vector to element face. The present work uses the numerical flux from the hybridized discontinuous Galerkin method developed by Samii *et al.* in (samii_etal_2019). Therefore, the numerical flux is defined through , an approximation to over the mesh skeleton called the numerical trace, as in (samii_etal_2019)

(18) |

where , and is a stabilization parameter motivated by the local Lax-Friedrichs numerical flux:

(19) |

In this description of the stabilization parameter,

is the maximum eigenvalue of the normal Jacobian matrix

:(20) |

The numerical trace must be such that the numerical flux is conserved across all internal edges in the mesh skeleton, and boundary conditions are satisfied at all boundary edges through the boundary operator defined according to an imposed boundary condition as in (samii_etal_2019):

(21) |

Eq.(17) and Eq.(21) form a system of equations that is used to solve for an approximate solution . For complete details of the formulation along with definitions for , see Samii *et al.* (samii_etal_2019).

In order to generate , a numerical solution operator for the dispersive correction part of the Green-Naghdi equations, Eq.(5) is written as a system of first order equations using the definition for operator (samii_and_dawson_2018):

(22) |

where . A variational formulation for Eq.(22) forms a global system of equations that would benefit from a dimensional reduction. Therefore, the hybridized discontinuous Galerkin method developed by Samii and Dawson in (samii_and_dawson_2018) is employed to treat numerically Eq.(22). According to (samii_and_dawson_2018), an approximate solution and are sought such that

(23) |

for all , where , and the numerical flux is defined as

(24) |

with a scalar constant used as the stabilization parameter. The numerical flux is weakly conserved and the imposed boundary conditions, defined through the boundary operator , are weakly satisfied as in (samii_and_dawson_2018):

(25) |

Eq.(23) is a series of local systems which forms block diagonal matrices that can be used to perform efficient static condensation of Eq.(25). This will form a global system of equations with its dimension equal to the dimension of . The system is solved to obtain that is subsequently substituted back into Eq.(23) to recover . The result is then used in the dispersive correction portion of the Green-Naghdi equations to seek an approximate solution that satisfies the variational formulation

(26) |

where . High order derivatives of , present in , are computed weakly using a discontinuous Galerkin method with centered numerical fluxes. See (samii_and_dawson_2018) for complete details of the presented formulation along with definitions of the boundary operators .

As a scalar conservation law, the Exner equation can be efficiently discretized using a discontinuous Galerkin method. To this end, an approximate solution is sought such that

(27) |

where a simple upwinding scheme is employed for the numerical flux since the sediment flux is not an explicit function of and the normal Jacobian matrix cannot be formed. Assuming that the sediment transport is always in the flow direction, the numerical flux is defined as (mirabito_etal_2011):

(28) |

where is the Roe-averaged velocity defined as

(29) |

In this description, superscript denotes a variable value at when approaching from the interior of an element , and when approaching from the exterior.

### 3.3 Coupled model

In the coupled model method, the Green-Naghdi and Exner equations, Eq.(2) and Eq.(8), are fully coupled and solved simultaneously. The Strang operator splitting technique is used also for the coupled model and the numerical solution operator for the dispersive correction part, , is as in the decoupled model. However, the operator has to be modified since it now needs to provide a numerical solution to the coupled system of the nonlinear shallow water and Exner equations and not only to the nonlinear shallow water equations as in the decoupled model.

The discontinuous Galerkin method developed for hyperbolic nonconservative partial differential equations by Rhebergen *et al.* (rhebergen_2008) is used, in the form presented by Mirabito *et al.* (mirabito_etal_2011), for the model that couples the nonlinear shallow water and Exner equations. In this method, the numerical scheme for the Exner equation is as in the decoupled model but the numerical scheme for the nonlinear shallow water equations requires corrections due to the nonconservative term present in the source term . Defining

, introducing a third order tensor

such that , and setting , we require that an approximate solution to the nonlinear shallow water equations satisfies the variational formulation (mirabito_etal_2011)(30) | ||||

where , is a Lipschitz continuous path from to such that and , and where with the superscripts and corresponding to elements and such that . The choice of the form for the path has minor effect on numerical solutions (rhebergen_2008); therefore, a simple linear path has been chosen for this numerical formulation. Subsequently, the integral in the nonconservative term may be evaluated as (mirabito_etal_2011)

(31) | ||||

It is worth noting that is single valued over the edges of the mesh skeleton and does not depend on the way the elements and are chosen for the edge . The numerical flux for this numerical scheme is defined as (rhebergen_2008)

(32) |

where the truncated characteristic speeds and are

(33a) | ||||

(33b) |

and the Harten–Lax–van Leer flux is (harten_etal_1983)

(34) |

Finally, the modified numerical solution operator seeks an approximate solution such that (mirabito_etal_2011)

(35) | ||||

### 3.4 Wetting-drying, wave breaking, and slope limiting

In the developed hydro-morphodynamic model, the water depth is assumed to be bounded from below by a positive value in the Green-Naghdi equations. This assumption must be ensured with a wetting-drying algorithm that preserves the positivity criterion for the water depth. Since the water depth is updated in the numerical solution operator only, the algorithm shall be executed in conjunction with the operator . In the presented work, the wetting-drying algorithm developed by Bunya *et al.* in (bunya_etal_2009) for discontinuous solutions to the nonlinear shallow water equations is adopted. Among the main features of the algorithm are: (1) the water depth is never allowed to drop below a specified minimum water depth , (2) the elements of the mesh used for numerical simulations are defined as ”wet” or ”dry” according to a classification algorithm, (3) the water mass is allowed to transfer from ”wet” to ”dry” elements only; otherwise, interfaces between ”wet” and ”dry” elements are treated as a reflecting boundary. For the dispersive correction and Exner equations a positivity preserving wetting-drying algorithm is not required, and the wetting-drying fronts are modeled as reflecting boundaries.

Although the Green-Naghdi equations have the ability to capture dispersive properties of water waves, the equations do not accurately resolve wave breaking phenomena in surf zones (bonneton_etal_2011). A more suitable depth-averaged hydrodynamic model capable of capturing wave breaking phenomena is formed by the nonlinear shallow water equations (bonneton_etal_2011). The use of the Strang operator splitting technique for the numerical treatment of the presented model provides an opportunity to switch between the Green-Naghdi and nonlinear shallow water equations in areas where one model is deemed to be more accurate than the other. In the developed splitting technique, it is possible to switch to the nonlinear shallow water equations by setting in regions where the Green-Naghdi equations can no longer provide an adequate approximation, e.g. in wave breaking regions. Therefore, a wave breaking detection criterion should be considered. To this end, the wave breaking criterion adopted by Duran and Marche in (duran_and_marche_2017) from the discontinuity detection criterion of Krivodonova *et al.* (krivodonova_etal_2004) is incorporated into the numerical model. The criterion states that wave breaking occurs over an element if the parameter (duran_and_marche_2017)

(36) |

is greater than a specified threshold that is typically . In this description of the parameter , is the element diameter, are the inflow faces of the element where , and is the total length of the inflow faces.

In applications of discontinuous Galerkin methods for the nonlinear shallow water equations, a slope limiter may be required in order to remove oscillations at sharp discontinuities in numerical solutions and preserve numerical stability. In particular, the wave breaking phenomena present themselves as sharp discontinuities in the numerical solutions. Therefore, the Cockburn-Shu limiter (cockburn_shu_2001) is incorporated into the numerical model and applied in conjunction with the operator . Changes in bed elevation may also form shocks that require a limiting procedure to avoid spurious oscillations in numerical solutions; thus, the Xu et al. limiter (xu_etal_2009) is integrated into the model to perform slope limiting in the Exner equation. The details of the limiters are not presented here, but readers are encouraged to consult the original sources.

## 4 Numerical experiment and discussion

The developed numerical model has been implemented in a software framework written in C++ programming language with the use of open source scientific computing libraries, such as Eigen

(eigen), Blaze (blaze), and PETSc (petsc). The software has been parallelized for shared and distributed memory systems with the use of a hybrid OpenMP+MPI programming, and HPX (hpx). Performance comparison between the hybrid programming and HPX has been performed by Bremer*et al.*in (bremer_etal_2019). The software has the capacity to simulate water waves using the discontinuous Galerkin finite element discretizations of the nonlinear shallow water and Green-Naghdi equations developed in (kubatko_etal_2006_nswe; samii_and_dawson_2018; samii_etal_2019), and it has been extended with the developed coupled and decoupled numerical models to allow for the possibility to simulate hydro-morphodynamic processes in coastal regions under the action of highly dispersive water waves.

The model has been validated against the experiment conducted by Sumer *et al.* (sumer_etal_2011) to measure flow and bed morphology induced by a solitary wave over a sloping beach. In the experiment four solitary waves have been run over a sloping beach inclined at rate; and, subsequently, a number of measurements have been performed, such as extents of sediment erosion and deposition over the sloping beach, and the free surface elevation at nine measuring stations. The measuring stations are located at the toe of the sloping beach, and at eight sections located 4.63, 4.69, 4.87, 5.11, 5.35, 5.59, 5.65, and 5.85 meters from the toe.

The choice of this experiment for the validation of the model has been motivated by the following reasons: (1) dispersive wave effects are prevalent in this experiment, and the Green-Naghdi equations should be used to resolve accurately the water wave dynamics, (2) in this experiment the solitary waves have sufficiently high amplitude to experience wave breaking; therefore, a wave breaking detection is required to switch to the nonlinear shallow water equations in surf zones, (3) in this experiment the sloping beach undergoes substantial sediment erosion and deposition that affect the bed surface elevation of the beach. Thus, performing numerical simulations of this experiment and comparing the results to the experimental ones have the potential to showcase all key features of the presented numerical model, such as the ability of the Green-Naghdi equations to simulate accurately water motion and capture dispersive wave effects, capacity of the numerical model to detect wave breaking regions and switch to the nonlinear shallow water equations in such regions, and the facility of the model to estimate sediment erosion and deposition due to bed-load transport.

To carry out the numerical simulations, a problem domain is partitioned into a finite element mesh comprised of square cells containing 2 triangular elements. The Dubiner polynomials of order are used for the approximating spaces (dubiner_1991). The sloping beach toe is located at , and all boundaries of the mesh are specified as reflecting boundaries. A two stage second-order Runge-Kutta method is used to perform time integration with a time step . The initial conditions for solitary waves in this experiment are characterized by equations (samii_and_dawson_2018)

(37) |

where the reference water level , the solitary wave height , the initial wave position , and

(38) |

Finally, for these simulations the bottom friction force is introduced into the numerical model through the source term by setting

(39) |

where the Chezy friction coefficient .

The simulations have first been performed over a rigid bed to validate the hydrodynamic model. Using separately the Green-Naghdi and nonlinear shallow water equations, the simulations have been run for 20 seconds which is a sufficient time for solitary waves to run up and run down along the sloping beach in this experiment. Fig.2 presents the free surface elevations obtained at the measuring stations from the experiment by Sumer *et al.* (sumer_etal_2011). As expected, in terms of accuracy the Green-Naghdi equations substantially outperform the nonlinear shallow water equations in the run up stage at the measuring stations located offshore. It is also evident that solitary waves break too early in the nonlinear shallow water equations simulations. In fact, the experimental results suggest that wave breaking occurs somewhere between the sections 3 and 5 which is accurately captured by the Green-Naghdi equations. However, neither model is able to accurately capture the water motion in the swash zone as evidenced by the free surface elevation measurements at the onshore section 8. We believe that these inaccuracies are due to the nontrivial physics that govern the water motion in swash zones, and to the limitations of the wetting-drying algorithm used in the simulations. Subsequently, the models are unable to capture correctly the water motion during the run down stage of the simulations. Nonetheless, the results are deemed satisfactory given the complexity of the physical processes occurring in flows induced by solitary waves over a sloping beach.

For the erodible bed simulations, the Grass model (grass_1981) in Eq.(10) for the sediment flux has been used with . The simulations with both coupled and decoupled model have been performed for 2 minutes and 30 seconds which is a sufficient time for water to substantially settle during the simulations. In each simulation four solitary waves have been run over the sloping beach. In this experiment the time scales in the hydrodynamic and morphodynamic models are comparable (li_etal_2019). Therefore, in the decoupled model the bed surface update has to be performed every time step of the hydrodynamic model. The bed surface erosion and deposition results obtained from the simulation runs are presented in Fig.3 for the decoupled model and in Fig.4 for the coupled model. As expected, the bed surface evolution in the offshore area is accurately estimated by both models since the hydrodynamic part captures the water motion in that area with sufficient accuracy. On the other hand, in the onshore area the models capture sediment erosion and deposition less accurately due, in part, to low accuracy of the hydrodynamic model in the swash zone. Overall, the results are considered satisfactory and indicate a promise for further development of the presented hydro-morphodynamic model, e.g. towards the extension of the model with suspended-load transport. Moreover, the decoupled model performed well relative to the coupled model and can provide a viable alternative to the coupled model, in particular, in cases where the time scales in the hydrodynamic part are shorter than in the morphodynamic part.

## 5 Conclusions

In this paper a hydro-morphodynamic model that couples a depth-averaged dispersive water wave model, the Green-Naghdi equations, with the Exner equation has been introduced. Although there are numerous works that couple the nonlinear shallow water equations with the Exner equation, to the best of authors’ knowledge, the coupling of the sediment continuity model with the Green-Naghdi equations has not been attempted before this work. The presented model is well suited for studying the bed surface evolution under bed-load transport in areas where dispersive wave effects are prevalent and should thus be included in the hydrodynamic model.

Numerical methods that utilise discontinuous Galerkin finite element methods have been presented for the hydro-morphodynamic model. The Strang operator splitting technique has been employed to single out the dispersive part of the Green-Naghdi equations for separate treatment. The resulting numerical models are augmented with wetting-drying, breaking wave detection, and slope limiting features. The numerical models have been used to simulate flow and sediment transport induced by solitary waves over a sloping beach. Comparing the numerical results with the experimental results collected by Sumer *et al.* (sumer_etal_2011), the numerical experiments have demonstrated that the presented model is capable of modeling water waves and sediment transport with a satisfactory accuracy in areas where the wave dispersion effects prevail.

The presented hydro-morphodynamic model shows a lot of promise. The depth-averaged dispersive wave model is not only able to accurately capture the water motion but also reduces the computational effort required to perform the simulations. The simulations carried out in this work have taken only a few minutes to run. Modeling hydro-morphodynamic processes in the same experimental setup using the Reynolds-averaged Navier-Stokes equations coupled with a sediment transport and morphodynamic model required simulations that took 3 days to run (li_etal_2019). Therefore, the presented model has a great potential to be used in simulations of hydro-morphodynamic processes caused by dispersive waves in large coastal areas.

## 6 Acknowledgments

This work has been supported by funding from the National Science Foundation Grant ACI-1339801, and the Portuguese government through Fundação para a Ciência e a Tecnologia (FCT), I.P., under the project DGCOAST (UTAP-EXPL/MAT/0017/2017). Authors would like to acknowledge the support of the Texas Advanced Computing Center through the allocation TG-DMS080016N used in the parallel computations of this work.