Macro-micro decomposition for consistent and conservative model order reduction of hyperbolic shallow water moment equations: A study using POD-Galerkin and dynamical low rank approximation

Geophysical flow simulations using hyperbolic shallow water moment equations require an efficient discretization of a potentially large system of PDEs, the so-called moment system. This calls for tailored model order reduction techniques that allow for efficient and accurate simulations while guaranteeing physical properties like mass conservation. In this paper, we develop the first model reduction for the hyperbolic shallow water moment equations and achieve mass conservation. This is accomplished using a macro-micro decomposition of the model into a macroscopic (conservative) part and a microscopic (non-conservative) part with subsequent model reduction using either POD-Galerkin or dynamical low-rank approximation only on the microscopic (non-conservative) part. Numerical experiments showcase the performance of the new model reduction methods including high accuracy and fast computation times together with guaranteed conservation and consistency properties.


Introduction
Accurate simulation of free-surface flows is necessary for prediction of natural hazards like floods, landslides, tsunami waves, and weather forecasting [14,15].However, the solution of the full incompressible Navier-Stokes equations is often too costly and simplified models like the shallow water equations (SWE) often yield inaccurate results, since the SWE model assumes a constant velocity over the vertical axis.The recently derived hyperbolic shallow water moment equations (HSWME) overcome this problem by allowing for polynomial velocity profiles [43].The model is based on a hyperbolic regularization of [46] using techniques from kinetic theory [24,44].The increased accuracy of the model was shown in numerical simulations, including sediment transport [26].However, the HSWME also lead to a higher computational cost and memory footprint due to additional nonlinear equations for the expansion coefficients of the polynomial velocity profile.The goal of this work is to perform an additional model reduction to improve efficiency of the model while at the same time preserving important properties like conservation of mass.
To reduce computational costs during the simulation of the HSWME we propose to project the system on a lower dimensional manifold using Galerkin methods.In this context we apply two model reduction techniques.First, the classical offline-online procedure (see for a review [48]), where important modes are calculated using a proper orthogonal decomposition (POD) in a so-called offline phase.This allows to conduct further online computations with the generated basis at significantly reduced computational costs.Second, the online adaptive basis method, dynamical low-rank approximation (DLRA) [39], that along with other methods like the adaptive basis and adaptive sampling discrete empirical interpolation method (AADEIM) [59], adapts the dominant modes by online updates to the solution locally in time.DLRA therefore does not require an offline computation since training is shifted to the online phase.
Most commonly in model order reduction (MOR) a reduced basis is computed with the help of the POD, that was first introduced in [55] in the context of fluid dynamics.In combination with Galerkin projection methods, POD-Galerkin has been successfully applied in fluid dynamics [51], including the SWE [53,52,71,64,65].Since POD-Galerkin does not preserve conservation laws, one possibility is to introduce Hamilton formulations that ensure structure-preservation when constructing reduced order models (ROM) of SWE with the help of POD [35,36].In this paper we consider the HSWME that introduce N additional equations to the 1D-SWE to account for non-constant velocity profiles along the vertical axis.With the aim to reduce the HSWME to a system yielding a similar complexity as the SWE we then apply POD-Galerkin as the first approach.In contrast to existing SWE-ROM approaches [53,52,71,64,65,35,36], that try to reduce the high dimensional state-space by Galerkin projections onto the spatial modes, we utilize POD modes that reduce the system only in the vertical direction.Here, the idea is to replace the N + 2 dimensional HSWME system by a r + 2 (where r ≪ N ) dimensional system.This is achieved through a linear combination of global ansatz functions that are determined by POD.Similar ideas have been first introduced by [67,68,69] and later in fluid dynamics to study the pipe flow where a coarse model featuring the dominant dynamics in the flow direction is enriched by a fine model that accounts for the additional dynamics in the transverse direction with the help of a modal expansion [60,62,61,56].This so called hierarchical model order reduction was further expanded for non-linear PDEs in [63].
Dynamical low-rank approximation [39] for matrix differential equations approximates the solution by a low-rank matrix decomposition and derives evolution equations to update the factors of the solution in time.These evolution equations are determined by minimizing the defect while restricting the evolution of the solution to the manifold of low-rank matrices.Stable time integrators for the resulting DLRA system, which are robust irrespective of the curvature of the low-rank manifold [38], are the matrix projector-splitting integrator [54] as well as the "unconventional" basis update & Galerkin step (BUG) integrator of [11].Here, we use the BUG integrator which only evolves the solution forward in time, thereby facilitating the construction of stable spatial discretizations [49].Moreover, the BUG integrator enables a straightforward basis augmentation step [9] which simplifies the construction of rank adaptive methods [9,50,31] and allows for conservation properties [9,22].
Comparing computational results for POD-Galerkin and DLRA, we note that the generation of basis functions during the online computation makes dynamical low-rank approximation significantly more costly than the online computation of POD-Galerkin.However, since DLRA does not require an offline computation in which the full model needs to be evaluated repeatedly, it offers three main advantages over POD-Galerkin: First, when the full order model is expensive to compute or cannot be stored in memory at a desirable accuracy, the offline phase of POD-Galerkin can become unfeasible to compute.Second, DLRA allows a straightforward application in settings that require a retraining of the basis with POD-Galerkin.Third, since DLRA allows to adapt the basis during the online phase it is independent of the choice of snapshots and it is beneficial in settings where a precomputed linear basis might not describe the dynamics sufficiently.In particular DLRA is superior to POD-Galerkin for transport dominated systems, where the Kolmogorov n-width decays slowly [58,27].
While model order reduction techniques promise a reduction of cost and memory, they can struggle to preserve important quantities of the full system such as boundary conditions or mass conservation.To guarantee preserving these properties, we propose to decompose the HSWME solution into two parts: (1) macroscopic variables: water height and momentum and (2) microscopic variables: higher-order moments of the velocity profile.Separate evolution equations for the macro and micro solution parts are derived and the model order reduction strategy is only applied to the computational and memory intensive microscopic part (2).Thereby, we facilitate imposing boundary conditions for the macroscopic water height and momentum while guaranteeing mass conservation.In addition, we achieve consistency with the underlying SWE model in case of vanishing microscopic structures.
The construction of efficient numerical method for dynamical low-rank approximation requires a carefully chosen formulation of the low-rank evolution equations and the construction of adequate numerical discretizations.Therefore, besides the introduction of a macro-micro decomposition to preserve solution invariants in model order reduction, main novelties of this work are: • The derivation of efficient low-rank evolution equations for HSWME.We derive evolution equations of low-rank factors of the HSWME which do not require the computation of full-rank matrices.Thereby, the memory footprint and computational costs to evolve the dynamical low-rank approximation are significantly reduced.
• The construction of an efficient numerical scheme.A classical implicit time discretization of friction terms yields prohibitive costs when evolving the low-rank approximation.We propose a splitting step to eliminate computationally costly terms in the numerical scheme, allowing for a significant reduction in runtime.
The remaining parts of this paper are structured as follows: We present the main concept of the macro-micro decomposition as well as the correction terms for model order reduction in Section 2. The mathematical model HSWME is briefly explained in Section 3. Section 4 covers the numerical discretizations including the implementation of the macro-micro decomposition for the new model.The conservative POD-Galerkin method is adopted for the considered model in Section 5, followed by the dynamical low-rank approximation in Section 6. Numerical results showing the performance of both methods are presented in Section 7 and the paper ends with a brief conclusion.

Concept
A core issue of model order reduction is the potential violation of crucial properties of the original problem.To preserve these properties we choose to first decompose the original dynamics into evolution equations describing the macroscopic (conservative) quantities on the one hand and microscopic (non-conservative) correction terms on the other hand.A reduced model is then derived only for the microscopic corrections terms, whereas the original macroscopic evolution equations remain unaltered.Thereby, basis functions which are important for overall solution properties remain unchanged, whereas the overall solution complexity that mainly arises from correction terms is reduced through model order reduction.The concept is visualized in Fig. 1.
Figure 1: Concept of the macro-micro decomposition MOR approach for HSWME.The HSWME is decomposed in macroscopic variables u and microscopic variables v.While the macroscopic variables are computed with a standard conservative scheme, the microscopic variables are treated by MOR techniques DLRA or POD-Galerkin.
We consider equations of the following type where A(q)∂ x q denotes the transport term and g ν (q) denotes the friction term, which depends on a friction parameter ν.Furthermore, the variable vector q = (u, v) contains both macroscopic variables u and microscopic variables v.
For the numerical solution, we perform two main steps: (I) A first order operator splitting to Eq. ( 1), in which we split the transport part from the right-hand side friction part.(II) A decomposition of the solution into macroscopic (u) and microscopic (v) variables.Note that the operator splitting only introduce first order errors in time and space, which are of the same order as the discretization errors in the sub-steps.The macro-micro variable decomposition does not introduce an error but a decoupled solution of the macro and micro step will do so.The procedure can be summarized by the following steps Step 1: transport Step 2: friction and the subsequent macro-micro decomposition Step 1a: macro transport Step 1b: micro transport Step 2a: macro friction Step 2b: micro friction Remark.Note that macro-micro decomposition has been used in a wide variety of applications in combination with model order reduction, see e.g.[20,19].In contrast to these works, we employ the macro-micro decomposition not to construct asymptotic-preserving schemes, but for two purposes: 1) Guarantee mass conservation and 2) ensure consistency with the SWE when the microscopic part tends to zero.Note that the latter also ensures that the model order reduction method does not reduce the approximation accuracy of water height and mean velocity.That is, these two quantities of interest can be represented accurately, independent of any chosen basis functions.The choice of macroscopic and microscopic variables for a general model is not predetermined and this could change based on the interpretation of the model or a spectral analysis potentially identifying a spectral gap, e.g., [1].
Remark.The notion of macro and micro variables originates from the limit of the SWME with respect to vanishing slip length.This models a perfect slip bottom such that only the interior friction remains and drives the "microscopic" higher-order moments to zero resulting in a constant velocity profile [34].The remaining water height h and velocity u m evolve according to the "macroscopic" shallow water equations.This justifies the decomposition into a macro part evolving on the equilibrium manifold and a remaining micro part.For more details, we refer to [34].
In the next sections, we will first introduce the full (coupled) model Eq.(1), then discuss the space time discretization of the operator splitting of Eq. ( 2) and Eq. ( 3), before deriving a detailed scheme for the macro-micro decomposition system Eq.(4) -Eq.(7), which can then be efficiently used for model reduction algorithms like POD-Galerkin and DLRA of the micro steps Eq. ( 5) and Eq. ( 7).

Models for Shallow Flows
In the following, we introduce two full models for shallow water flows: the shallow water moment equations and their hyperbolic formulation, the hyperbolic shallow water moment equations.
We note that an alternative approach to starting with these models is to directly apply the model reduction techniques to the underlying incompressible Navier-Stokes equations.Conceptually, this is possible.In the case of DLRA, the evolution equation for the L-step is then continuous in the velocity variable, and a moment method can be used to obtain a discretized set of equations.This is equivalent (when using the BUG integrator as we do in this work) to first performing the moment approximation and then applying DLRA.However, using the moment method first and the model reduction second is much more convenient in terms of notation and efficiency of presentation.It also allows for a very intuitive understanding of how to carry over the desirable properties of the moment model to the numerical solution.

Shallow water moment equations
Shallow water flows are often modeled using the standard SWE, which are hyperbolic balance laws for the macroscopic variables water height h and the mean velocity u m , derived from the underlying incompressible Navier-Stokes equations.For simplicity, we consider a one-dimensional vertical coordinate x ∈ R, a flat bottom topography, and a Newtonian fluid.The SWE model for h and u m then reads [46] with gravity constant g, slip length λ and kinematic viscosity ν.Note that the system is written in so-called convective variables (h, hu m ), which will also be the case for the more advanced shallow water Moment model below.The addition of a non-flat bottom topography term in Eq. ( 8) will pose no conceptual difficulties and is left for future work, see [26,42] for examples of the straightforward treatment.Further note that the concept can be readily extended to the 2D case.However, so far there is very little research on 2D models with vertically resolved velocity profiles.Even though the conceptual framework of the model is laid out in [46] and the hyperbolic regularisation below can be extended based on [43], to the best of our knowledge, no full 2D simulations exist.We therefore focus on the 1D case in the current work.Due to its simplicity with only the macroscopic mean velocity u m , the SWE model cannot represent vertical variations of the velocity, representing micro structure of the profile.The assumption of a constant velocity profile breaks down for many applications, especially when considering strong bottom friction [37].Even when starting with constant velocity profiles, bottom friction leads to a deceleration of the fluid close to the bottom giving rise to more complex velocity profiles.
The recently developed shallow water moment equations (SWME) [46] tackle this problem by introducing a polynomial expansion of the velocity profile u(t, x, z) depending on the vertical variable z ∈ [0, h] as follows where ϕ j : [0, 1] → R are the scaled Legendre polynomials of degree j defined by The basis functions ϕ j form an orthogonal basis, due to , where δ mn is the Kronecker delta.
The corresponding expansion coefficients α j for j = 1, 2, . . ., N of the polynomial expansion Eq. ( 9) are also called moments in analogy to moment models from kinetic theory [66].The moments represent micro variables that augment the macroscopic variables h and u m .
The expansion Eq. ( 9) then allows for the representation of complex velocity profiles, e.g., linear, quadratic, ..., up to the maximal polynomial degree N ∈ N.While a larger maximal degree N leads to potentially higher accuracy of the velocity profile, more moments are used in the representation which leads to higher computational costs.
Evolution equations of the moments in time and space can be derived by a projection of the underlying incompressible Navier-Stokes equations onto the basis functions Eq. (10).In combination with the conservation of mass, the resulting system of equations is called the shallow water moment equations.More details on the derivation of the SWME can be found in [46].
In this paper, we do not use the SWME model in the form derived in [46] due to its lack of hyperbolicity, which was attributed to stability issues in [43].We therefore do not show the explicit form of the matrix A(q) and source term g ν (q) of the SWME here, but directly introduce the hyperbolic regularization derived in [43].

Hyperbolic shallow water moment equations
Based on results from kinetic theory [24,44,6], the HSWME were first derived in [43] and have been extended and analyzed in [42,34].As a key property of transport dominated problems, hyperbolicity is the property of the transport matrix to have real eigenvalues corresponding to waves with real propagation speeds.It was shown that the HSWME yield accurate results while preserving hyperbolicity of the model equations [43].Recently, the model has been applied to sediment transport [26] and analysis of steady states and equilibrium stability has been performed [34,42].
In this paper, we use the HSWME model in its standard form written in terms of the convective variable vector q = [h, hu m , hα 1 , . . ., hα N ] T ∈ R N +2 as with hyperbolic transport matrix A(q) ∈ R (N +2)×(N +2) given by Note that the transport matrix A is a function that only depends on macroscopic h, u m , plus α 1 , whereas it does not depend on the higher microscopic moments α i for i = 2, . . ., N .This will be the first important ingredient allowing for efficient model reduction via POD-Galerkin in Section 5 and DLRA in Section 6 later.While the HSWME model appears simpler than SWME, it still retains the main nonlinearity of the SWME and is able to reproduce the complex nonlinear flow patterns of shallow flows, as demonstrated in numerous recent papers [26,34,42,1].
The source term g ν (q) ∈ R N +2 of Eq. ( 11) is given by g ν (q) = [0, g ν 0 , . . ., g ν N ] T and reads with slip length λ, viscosity ν, and constants a i,j given by Note that the first entry of the right-hand side friction term g ν (q) is zero, leading to the conservation of mass, which is the integral of the water height h.More importantly, the other entries given by g ν i (q) only depend non-linearly on the water height h, while depending linearly on the remaining u m and α i .This will be another important ingredient for efficient model reduction to be exploited later.

Space-time discretization with operator splitting
The right-hand side friction term of Eq. ( 11) can be stiff for small λ or small h.This prohibits the use of standard explicit schemes to solve the coupled model Eq.(11).While stable explicit schemes like projective integration exist [1], they typically require more time steps and some parameter tuning, which would be impractical for a robust model order reduction later.The left-hand side transport part is naturally discretised best with an explicit scheme using a standard CFL condition.We therefore split the full model Eq. ( 11) into a transport step and a friction step.We then treat the space-time discretization of both terms separately.As explained in Section 2, for a single time step of a numerical solution we apply a first order operator splitting to Eq. ( 11) in the same fashion as [34], in which we split the transport part Eq. ( 2) from the right-hand side friction part Eq. (3) as Step 1: transport where the solution of the transport and friction step will be considered separately in the following two subsections.

A new macro-micro decomposition scheme for the transport step
A straightforward discretization of the transport step Eq. ( 2) via a standard path-conservative numerical schemes [42] is sufficient to solve the transport step of the full order model.However, preparing for the model reduction in Section 5 and Section 6, we here describe a new version exploiting the structure of the model using a macro-micro decomposition of the variables.This means additionally decomposing the vector of variables in the same way as they will be treated by the POD-Galerkin and DLRA methods in Section 5 and Section 6, respectively.A standard numerical scheme as written in Appendix A, see Eq. ( 24), considers an update of the full state vector.However, this does not allow to completely leverage the structure of the underlying model equations during the model reduction procedure later.To that end, we exploit the structure of the HSWME model Eq.(11), by decomposing Q into two parts: (1) the first two (macroscopic) variables for the water height h and the mean velocity u m (called U) on the one hand and (2) the last N (microscopic) variables for the moments α i on the other hand (called V).This can be written as In the same fashion, we decompose the transport matrix A(q) into four blocks corresponding to the first two equations and the last N equations and variables, respectively with blocks where all other entries are zeros.We thus decompose the system into a set of macroscopic variables U and another set of microscopic variables V .While the macroscopic variables are not reduced to allow structure preservation, the microscopic variables are used to generate an efficient reduced order model.
The space-time discretization of the transport part towards a solution scheme is described in the Appendix A for brevity.The derivation of reduced model by both POD-Galerkin as well as DLRA rely on the definition of this discretized scheme.
After the solution of the transport step has been computed the scheme continues with the friction step, which is potentially stiff and requires a different, implicit scheme.

Friction step
To solve the space-homogeneous friction step Eq. ( 3), an implicit scheme is necessary due to potential stiffness originating from small values for λ or h.We will use the scheme from [34], but adopt it to our decomposition of variables Q into U and V as explained in this section.
The macro variables u = (h, hu m ) after the friction step are obtained using the definition of the source in Eq. ( 12) and the two observations: (1) The height h remains constant in time during the friction step.(2) The micro moments α i remain constant during the macro update.This leads to the updated values (h n+1 , h n+1 u n+1 m ).For the update of the remaining micro coefficients V during the micro friction step Eq. ( 7), we make use of the fact that h and u m are now constant and V only occurs linearly.According to [34], Eq. ( 3) for the last N moments can then be written as where ).The micro friction step Eq. ( 15) is then solved implicitly by a backward Euler method to overcome stability issues from potential stiffness.Note again that the water height h j is constant during the whole friction step and u n+1 m,j is given by the solution of the macro friction step Eq. ( 6).Applying the backward Euler method to Eq. ( 15) yields We then define the matrix to arrive at the time update We note that D −1 j can be precomputed efficiently, as it does not depend on u m and α i .Thus, the updated micro variables in Eq. ( 17) can be computed efficiently without inverting a matrix during the online computation.

Macro-micro decomposition for conservative POD-Galerkin
The following section addresses the POD-Galerkin reduction of the HSWME model from Section 3.2.As already pointed out in Section 2 the presented approach applies the reduction only to the microscopic higher-order moments v = [hα 1 , . . ., hα N ] ⊤ ∈ R N .Furthermore, in contrast to the conventional POD-Galerkin approach that reduces the full state-space of the discretized PDE, we only reduce the dimensions of the PDEs' state-components.This leads to a decoupling between the time-space dynamics and the correlation between the components and avoids the separation of spatial-temporal dynamics that can lead to slow decaying approximation errors, known as the Kolmogorov N -width problem [58,27].For a typical dam break scenario we showcase this slow decay in Fig. 3 of our numerical Section 7. To obtain a reduced model that allows rapid yet accurate predictions over a range of different parameters, one often uses a two phase offline-online procedure, which is explained in the following.

Offline phase
In the offline phase the reduced basis and the operators of the reduced model space are precomputed.Note that in contrast to the DLRA approach from Section 6 one is willing to forego large up-front offline costs in favor of a more efficient reduced model for the online phase.The reduced basis is formed by collecting snapshots of the solution of Eq. ( 11) at N t different time/parameter instances V = {v(x, t 1 ), . . ., v(x, t Nt−1 )} .
The POD approximates the microscopic higher moment vector using an orthonormal basis where W = [w 1 , . . ., w r ] ∈ R N ×r collects the basis vectors and v = [ α1 , . . ., αr ] ⊤ ∈ R r denotes the reduced space-time coefficients.In a fully discrete setting this basis is computed by a truncated singular value decomposition (SVD) of the snapshot matrix where the V n are defined in Eq. ( 13).The truncated SVD of V POD yields Here, Σ = diag (σ 1 , . . ., σ r ), r ≪ M , is a diagonal matrix containing the largest r singular values A common choice for r is to truncate after a certain energy percentage (e.g.E cum ≥ 95%) is reached in the reduced system compared to the full system: Calculating an SVD can be challenging, particularly when dealing with higher-dimensional space, due to the potentially prohibitive size of the snapshot matrix.This can be circumvented by either using randomized or wavelet techniques to calculate the POD-modes [47,30] more efficiently or applying POD-greedy sampling methods [29,70] that reduce the number of snapshots needed to form a POD basis.
In the context of Galerkin projections the column space of W yields the trial space

Online Phase
As explained in Section 2, the online phase evolves the dynamics in two steps.In the first step we evolve the first two macro variables U without any model order reduction and in the second step we use the evolved U as an input of our reduced system of the microscopic higher order moments V.The reduced system is gathered by projecting Eqs. ( 16) and ( 25) onto the test space which is the reduced subspace spanned by the columns of W.

Transport step
For the transport part in Eq. ( 25) we define the reduced transport term as Evaluating this at the jth cell we obtain: where we define the reduced non-linear operators Âj,vu := W ⊤ A j,vu and Âj,vv , where: with identity matrix I r ∈ R r×r .Note, since Â is precomputed in the offline phase, the non-linear term Eq. ( 19) only has to be evaluated for r instead of N components, which generates the speedup.To further simplify the non-linear terms POD is often used in combination with sparse sampling methods that sample the nonlinear terms at a few components to approximate them in a low-dimensional space.Examples are the discrete empirical interpolation method (DEIM) [13], that relates back to the empirical interpolation method (EIM) [3,28], the gappy POD [23,5,2], or the energy-conserving sampling and weighting (ECSW) method [25].However, we refrain from using them here for a direct comparison with DLRA, that does not make use of sparse sampling methods, yet.

Friction step
Similar to the transport term we evaluate the micro friction term Eq. ( 7) at the jth component projected onto the test space.We define the components of the reduced source term Ĝ( V The complexity of the source term computation is reduced, since only r components have to be evaluated.

Dynamical low-rank approximation
As an alternative to POD-Galerkin, we propose to evolve the microscopic higher-order moments with DLRA introduced in [39].This method is data-driven in the sense that it closes the moment equations without assuming physical properties, based on the real-time solution data.

Macro-micro decomposition for dynamical low-rank approximation
The core idea of DLRA is to evolve the solution on a low-rank manifold.That is, DLRA represents the solution as a low-rank factorization and provides evolution equations for the individual factors.Therefore, DLRA can be interpreted as a Galerkin method which updates not only the expansion coefficients but also the basis functions in time.To preserve the structure of water height and momentum, we apply DLRA to the microscopic correction terms v i := hα i only.Collecting the spatially discretized correction terms in a matrix V(t) ∈ R Nx×N , where v ji = h(t, x j )α i (t, x j ), we define a low-rank approximation as V(t) = X(t)S(t)W(t) ⊤ , where X ∈ R Nx×r and W ∈ R N ×r can be interpreted as the collection of r basis vectors in space and moment order with a corresponding coefficient matrix S ∈ R r×r .
To benefit from this representation, we want to ensure that the method works on these factors only and never needs to compute and store the full solution V. To preserve the low-rank structure of the solution, we therefore force the solution at all times t to remain in the manifold of rank r matrices, which we call M r .This can be ensured when the time derivative of V lies in the tangent space of rank r matrices, i. e., the solution when advancing in time does not leave the manifold M r .We denote the tangent space at V(t) as T V(t) M r .Then, the time evolution equations for the basis vectors and coefficients must satisfy The first condition conserves the representation V(t) = X(t)S(t)W(t) ⊤ and the second condition minimizes the residual, which is essentially a Galerkin condition.Hence, an equivalent formulation of (20) is This represents a Galerkin method, which chooses test functions based on the solution data and the geometry of the manifold of rank r matrices.An admissible choice of δV when the solution reads V(t) = X(t)S(t)W(t) ⊤ is δV = X i W ⊤ j which indeed lies in the tangent space.Here, X i and W j denote the i th and j th columns of the basis matrices, respectively.This test function yields an evolution equation for the coefficient matrix Choosing test functions δV = X i and δV = W i , which again lie in the tangent space, results in evolution equations for the respective basis matrices Solving these evolution equations with classical time integration schemes is inefficient, since the coefficient matrices are often ill conditioned and the time step size is dictated by the smallest absolute eigenvalue of S. Two robust time integrators, which guarantee stability, are the projector-splitting integrator [54] as well as the BUG integrator [12].Their main strategy is to not evolve X and W in time directly, but to evolve a linear transformation K := XS and L := WS ⊤ and retrieve the basis matrices through a QR-decomposition.In this work, we focus on the BUG integrator, which consists of three update steps 1. K-step: Update X 0 to X 1 via .
The solution at the next time step is then given by V(t 1 ) = X 1 S 1 W 1,⊤ .Note that the first two equations which evolve the basis in time can be updated in parallel followed by a serial update of the coefficient vector.Effectively the size of the original matrix ODE is reduced from a N ×N x system to three smaller matrix ODEs of size: N x ×r (K-step), N ×r (L-step), and r×r (S-step).We wish to underline the requirements of a user-determined choice of the rank r, which needs to be determined such that the solution exhibits a sufficiently accurate approximation at minimal computational costs and memory requirements.Note however that the BUG integrator allows for an extension to rank adaptivity [9].Here, the basis after the K and L steps is augmented with the basis at time t 0 to enlarge the approximation space from r to 2r basis vectors.After a 2r×2r coefficient update in the S-step, the solution is truncated to a new rank r 1 ≤ 2r according to a user-determined tolerance parameter ϑ.While this parameter needs to be chosen before the simulation, it has a clear interpretation as the truncation error in every time step.For clarity of presentation, we focus on the fixed-rank BUG integrator in our derivations, the numerical results however include both fixed-rank and rank-adaptive computations.For more information on rank-adaptivity and further properties of the rank-adaptive integrator, we refer to [9].Further approaches for rank-adaptivity in dynamical low-rank approximation and model order reduction are, for example, [33,10,16,32,18].The classical dynamical low-rank approximation approach often does not preserve important physical properties.This stems from the fact that DLRA can remove basis functions which are needed for conservation.However, problem-dependent adaptations to the classical DLRA integrators can provide conservation properties.As an example, conservation of solution invariants up to a tolerance parameter can be achieved with a basis augmentation step [9].Moreover, [22] uses a basis augmentation as well as a reformulation of the K, L and S-step to preserve mass, momentum and energy in the Vlasov equations.
Our approach in this work ensures local mass conservation by decomposing the dynamics of the macroscopic (conserved) water height and momentum from the dynamics of the microscopic correction terms.This means that the conservation law structure of the water height equation is not altered by DLRA.

Evolution equations for the low-rank HSWME
In the following we derive efficient representations of the K, L, and S-steps for the micro transport step Eq. ( 5) and the micro friction step Eq. ( 7) of the HSWME model.Note that the macroscopic variables are not altered by the DLRA and computed during the macro transport and macro friction step as explained in Section 4. Note that the derivation of evolution equations is performed on the semi-discrete level, i.e., the system is already discretized in the vertical direction z by a moment approximation and the horizontal direction by a finite volume discretization.DLRA can, however, also be derived in the fully continuous formulation [21], leading to a continuous set of low-rank evolution equations.When using BUG integrators combined with an equivalent moment and finite volume discretization of the continuous low-rank evolution equations, these two approaches are equivalent, and we choose the discrete formulation for ease of presentation.

Transport step
For the micro transport part Eq. ( 5), the spatially discretized right-hand side Eq. ( 25) is given by (A j,vu (U j − U j−1 ) + A j,vv (V j − V j−1 )) .

K-step:
To derive the K-step, we denote the j-th row of K as K j (t) ∈ R r and represent the solution V at spatial cell j as V j (t) = W 0 K j (t).Multiplying the right-hand side F v with W 0,⊤ then yields where A j,vv := W 0,⊤ A j,vv W 0 .Note that since A j,vv ∈ R r×r and W 0,⊤ A j+1,vu ∈ R r×2 , the main memory requirements stem from storing K j at all spatial cells j, i. e., memory requirements are O(r • N x ) opposed to the original method's requirements of O(N x • N ).
An efficient computation of A j,vv precomputes A := W 0,⊤ AW 0 which requires O(r 2 • N ) operations.Using the matrix A ∈ R r×r to compute A j,vv via requires O(r 2 •N x ) operations.Hence, the computational costs for the K-step are C K ≲ r 2 •(N x +N ).

L-step:
To derive the L-step, we represent the solution V at spatial cell j as V j (t) = L(t)X 0 j .Moreover, we test with X 0 , meaning that we multiply the right-hand side with X 0 j and sum over j.Hence, using Einstein's sum convention and writing out A j,vv we have To simplify the structure of the resulting equations, we define and α 0 − as well as u 0 − accordingly.This results in Note that the highest memory requirements come from storing L ∈ R N ×r .Computing X 0 , α 0 ± and u 0 ± requires O(r 2 • N x ) operations and computing A • L requires O(N • r) operations, since A only has off-diagonal entries.The multiplication of L and X 0 , α 0 ± as well as u 0 ± requires O(r 2 • N ) operations, i. e., we have a computational cost for the L-step of C L ≲ O(r 2 • (N x + N )).

S-step:
To derive the S-step, we represent the microscopic solution V at spatial cell j as V j (t) = W 1 S(t) ⊤ X 1 j .Moreover, we test with X 1 and W 1 , i. e., we multiply the right-hand side with X 1 j and sum over j as well as multiply with W 1,⊤ .With the previous definitions we have The memory requirements are O(r 2 ) and following the discussion for the K and L-steps, the computational costs for the S-step again are C S ≲ O(r 2 • (N x + N )).

Friction step
For the friction step Eq. ( 7) we have according to Eq. ( 15) For simplicity of notation, we define the matrix ordinary differential equation (ODE) where ) will be the values obtained from a previous macro friction update Eq. ( 6) and assumed known during the micro step Eq. ( 7).Let us again derive K, L and S-steps and directly define a time discretization.

K-step:
Let us use the representation V = K(t)W 1,⊤ and test with W 1 .Then, the K-step equation reads Since friction terms are commonly stiff, we use an implicit Euler method to discretize in time.That is, for a fixed spatial cell j and defining where u n+1 m,j is the updated velocity from the macro friction step Eq. ( 6).We again retrieve the time updated spatial basis which we denote by X 2 by a QR decomposition of K n+1 .Hence, at N x spatial cells, we need to solve a linear system with r unknowns.Moreover, we have to compute flux matrices with r 2 entries which require N operations per entry.Thus, the computational costs are L-step: Use the representation V = X 1 L ⊤ and test the transposed of ( 21) with X 1 .Then, the L-step equation reads Let us define h −2 1 := X 1,⊤ h −2 X 1 and h −1 1 := X 1,⊤ h −1 X 1 and again use an implicit Euler time discretization.This gives The time updated coefficient basis which we denote by W 2 is then obtained by a QR decomposition of L n+1 .Hence, we need to solve an r • N system of linear equations which requires O(r S-step: Use the representation V = X 2 SW 2,⊤ and test ( 21) with X 2 and W 2 .Then, the S-step equation reads An implicit Euler time discretization gives with Hence, we need to solve an r 2 system of linear equations which requires O(r 6 ) operations.Moreover, computing flux matrices requires O((N x + N ) • r 2 ) operations, hence Remark.Note that the implicit L-step exhibits high computational costs of O(r 3 • N 3 ).Nevertheless, we are able to separate the spatial degrees of freedom from the number of moments, which are commonly much smaller than the number of spatial cells.Therefore, we expect a significant reduction of computational costs by DLRA, which, however, does not resemble its full potential.Note that DLRA will become substantially more efficient compared to the full-rank baseline for two-dimensional spatial domains since, in this case, the number of spatial cells drastically increases and thereby removes the computational bottleneck of matrix inversions performed in the L-step.In this case, the computational costs are distributed more evenly among the K, L, and S steps.

Numerical experiments
As numerical examples of the model reduction techniques developed in the previous sections, we consider three test cases: (1) a water column (also called dam break) test case similar to the test case in [43] (2) a smooth wave similar to the test case in [46] and (3) the square root velocity profile test case that mimics a realistic velocity profile.Additionally to POD-Galerkin and DLRA we also include computations of the HSWME with a reduced number of moments (rHSWME).In the following, r denotes the rank in the case of DLRA and POD-Galerkin and the number of Moments for the rHSWME.All numerical experiments can be reproduced with the openly available source code [45].The simulations have been performed on 11th Gen Intel(R) Core(TM) i7-11850H CPUs (8 processing units).

Water column or dam break test case
In the first test case, we investigate a water column or dam break scenario defined by the settings in Table 1.Initially, the water is at rest and the height profile is defined as h(x) = 0.3 + 0.35 • (tanh(x) − tanh(x − 0.2)), for x ∈ [−1, 1] leading to a water column within [0, 0.2] with slightly smoothed boundaries.Note, that we are using a smoothed water column, because of the oscillator behavior of a Lax-Friedrichs scheme with non-smooth initial conditions [4].The friction parameters are ν = 1.0 and λ = 0.5 which slow down velocities close to the bottom.
The numerical discretization is performed using N x = 2000 cells and a CFL number of 0.25 for time stepping within t ∈ [0, 0.2].For the full-order HSWME model, we consider N = 100 coefficients, leading to a large system of coupled PDEs.For the model reduction, we first consider a fixed number of r = 5 basis function, which largely reduces the complexity of the reduced POD-Galerkin and DLRA models.
The trial and test space for the POD-Galerkin approach is spanned by the POD-basis.It is setup offline from the concatenation of the snapshots for two different trajectories sampled with slip length ν ∈ {0.1, 10}.The snapshots are shown for the water height in Fig. 2. From the snapshot matrix V POD ∈ R (2NtNx)×N we can compute the POD basis friction coefficient [42] Table 1: Simulation setup for water column test case.
with help of the SVD Eq. ( 18) and evaluate the projected right-hand side efficiently for the testing slip length ν = 1.We highlight that our approach does not try to decouple space and time and we therefore avoid slowly decaying approximation errors of our dyadic decomposition, a known problem for transport-dominated fluid systems.The difference between the classical MOR approach, which suffers from this slow decay and the here presented approach is shown in Fig. 3.As our approach only tries to reduce the component space, it avoids slowly decaying approximation errors, thus yielding a rapid decay of the POD approximation errors.Note, that this procedure is nothing new but has been extensively used in hierarchical MOR reduction (see for example [63]).3: Relative approximation error i ∥v(x, t i ) − ṽ(x, t i )∥ L 2 of the snapshot data for space-moment ṽ(x, t) = k α k (t)w k (x) and moment only basis ṽ(x, t) = r k=1 αk (x, t)w k .The online error is given as ∥u(x, t end ) − ũ(x, t end )∥ L 2 , where u = (h, hu m ) basis functions are chosen.This emphasizes the good approximation quality of the reduced models despite a small number of basis functions.However, we note that for this specific test case the moment model with only 5 moments already does a good job as well.This is because the dependencies on the higher moments vanishes quickly and a small number of moments already gives an excellent approximation quality.In Fig. 5 the water velocity profiles are plotted each time for different positions: x = 0.05, 0.0, 0.15 close to the center of the domain in Fig. 5a and x = 0.65, 0.67 close to the shock wave in Fig. 5b.In Fig. 5a, we clearly see that also the velocity profiles of the full-order model and the reduced POD-Galerkin and DLRA models, as well as the reduced moment model agree almost perfectly at all three points.The SWE model on the other hand shows an overestimation of the average velocity at both positions.This is due to wrong propagation speeds of the SWE model [43] and clearly shows why this simple model is not useful in simulations of such complex cases.In Fig. 5b, the profiles close to the shock wave are plotted and the reduced models again agree almost perfectly with the full-order model.Again the SWE model predicts wrong average profiles, in this case much smaller than the full model and the reduced models.
For the settings from above, Fig. 6 shows the runtime comparison between the full-order HSWME, DLRA, SWE, and POD-Galerkin, where the POD-Galerkin runtime is divided into the offline precomputations and the online phase.The truncation ranks of DLRA r = 4 and POD-Galerkin r = 3 are tuned such that the relative L 2 errors of u at the final time are approximately the same (≈ 0.3%, compare Fig. 7).The DLRA method already reduces the runtime significantly.The SWE has the fastest runtime, but does not result in sufficient accuracy as seen in the previous figures.The POD-Galerkin method is equally as fast as the SWE during the online phase, but requires a relatively costly offline precomputation phase.It is important to highlight that the efficacy of a global reduced order model becomes particularly apparent in the realm of multi-query simulations.In this context, the upfront costs associated with the offline phase of the POD-Galerkin algorithm are offset by an exceptionally efficient online phase, which can be replicated across a considerable number of queries.
However, there might be other MOR applications in which a global ROM can not be set up, because of the tremendous amount of data that would be required in the offline stage.For example, in the realistic scenario of two-dimensional HSWME or higher-dimensional systems involving multiple parameters formulated in a tensor-valued ODE (discussed in [40]).This is where DLRA has a clear advantage, as it does not rely on an offline phase.When the offline phase becomes too expensive, a combination of both methods might be beneficial.In this approach, DLRA is used to establish the initial basis, which is then employed in the online phase of POD.
Increasing the rank of the reduced models leads to increasing accuracy at the expense of more runtime, as can be seen in Fig. 7, where we compare the online computation time of POD-Galerkin with DLRA and the reduced moment HSWME model.For the rank adaptive version of DLRA we sample the solution for 14 different ϑ values evenly on a log scale in the interval 10 −10 ≤ ϑ ≤ 10 0 .From the comparison we observe, that rHSWME is the fastest, which is not surprising since the higher moments vanish quickly, in a regime of high friction.POD-Galerkin seems to achieve similar errors for small number of modes, however, it is slower then the rHSWME with increasing accuracy.Both DLRA and its adaptive version are one order of magnitude slower than POD-Galerkin.Note that since we apply MOR only in the microscopic higher moments (v-components) the reduced DLRA model converges towards the SWE with r → 0. Therefore, we obtain small relative errors even with a small number of modes.
As a summary of this first test case, we see very good accuracy of the reduced models with small rank and significant speedup that can even be amplified by further reducing the rank.Unfortunately, these speedups do not pay off directly when comparing it to the reduced moments HSWME, as the higher moments vanish quickly, due to the choice of the test case.In the next section, we will see, that this is not generally the case.

Smooth wave
This second test case closely follows the general simulation test cases in [46] and its setup is given in Table 2.It describes a smooth wave given by the initial height function h(x) = 1 + exp(3 cos(π(x + 0.5)))/exp(4) travelling through a periodic domain The initial velocity profile is chosen as u(0, x, ζ) = 0.25 , leading to u m = 0.25, α 1 = −0.25, and α N = 0.25.This velocity profile happens to be outside of the hyperbolicity region of the standard SWME model in the test case in [43].This means that applying the SWME model can lead to instability issues.This is resolved by applying the HSWME model and demonstrates the utility of this guaranteed hyperbolic model necessary to result in a well-posed model.Given the substantial dependence of the initial profile on the last moment α N , it is anticipated that the rHSWME will be surpassed by DLRA and POD-Galerkin.Similar to the water column test case, we compute the POD-basis from the snapshots of two trajectories simulated with the full HSWME at ν ∈ {10, 1000}.The friction parameters to compare all methods are now chosen as ν = 10 and λ = 0.001, leading to large values of the bottom friction term.We first show the numerical results for the macroscopic water height h and momentum hu m in Fig. 8 for (1) the full-order HSWME model, (2) the lowest-order SWE model, (3) the reduced HSWME model, (4) the new macro-micro decomposition POD-Galerkin, and (5) the new macro-micro decomposition DLRA.As for the first test case, the simple SWE model fails at capturing the complex dynamics of the test case, which is especially apparent in the momentum hu m shown in Fig. 8b.Similarly, POD Galerkin and the rHSWME model struggle to capture an accurate description of the momentum.Here DLRA yields a very good match even with a small number of degrees of freedom.
In Fig. 9 the water velocity profiles for the smooth wave test case are plotted each time for different positions: x = 0.05, 0.0, 0.15 close to the center of the domain in Fig. 9a and x = 0.65, 0.67 further away from the center in Fig. 9b.In Fig. 9a, the velocity profiles of the full-order model agree well with the reduced POD-Galerkin and DLRA models at all three points, while there are some small differences in the maximum velocity value due to the complexity of the test case and the significant model reduction.In contrast the rHSWME completely overshoots in the maximum velocity as it misses the information of the highest moment.Similarly for the SWE model, that is not shown for conciseness.The velocity profiles at positions further away from the center in Fig. 9b yield a similar result.For this smooth wave test case, we also want to emphasize one main property of our newly developed macro-micro decomposition reduced models, which is guaranteed mass conservation by construction.In Fig. 10, the time evolutions of the total mass 1 −1 h dx, the total momentum 1 −1 hu m dx and the total fourth higher moment 1 −1 hα 4 dx are plotted in terms of relative deviation from the initial values.It is clearly seen that the total mass is constant as changes are within machine precision.The total momentum and total higher momentum are not conserved in agreement with the underlying PDE model, which includes friction terms for the corresponding equations leading to a loss of momentum due to bottom friction, for example.The mass conservation is achieved due to the macro-micro decomposition formulation which includes explicitly solving for the water height while applying the model reduction only to the remaining microscopic velocity profile coefficients.While the total momentum is not constant, its evolution is in very good agreement with the full HSWME model for both reduced models since also the equation for hu m is apart from the model reduction process.Some deviations are seen for the total fourth momentum, obviously originating from the model reduction process for that variable.
With Fig. 11 we want to emphasize that a naive application of model reduction techniques does not lead to conservation of mass.This is done by comparison of our new macro-micro decomposition conservative DLRA, where the evolution of the macroscopic water height h and momentum hu m is decoupled from the microscopic reduced coefficient system, with a naive (non-conservative) DLRA, where the complete system including water height h, momentum hu m and coefficients hα i is reduced as a whole.It is seen that only the macro-micro decomposition conservative DLRA method achieves mass conservation.
Next, we compare the runtime vs. speedup of the presented methods.For this study, the rank of DLRA and POD-Galerkin is gradually increased starting from r = 1 until r = 30 (not all points in between are included).As before the rank adaptive version of DLRA is sampled for 14 different ϑ values evenly on a log scale in the interval   10 −10 ≤ ϑ ≤ 10 0 .Comparing both MOR methods and the rHSWME in Fig. 12, we observe that rHSWME stagnates at about 2% error.This can be explained as after the linear moment the higher moments only increase the complexity of the calculations, but can not capture any sensible dynamics before the last moment is not included.DLRA and POD-Galerkin however do not a-priori impose a basis and therefore can include the information necessary to represent the dynamics of the highest moments.As expected, we observe that DLRA exhibits an increased runtime in comparison to the online phase of POD-Galerkin.In general, we see that the adaptive DLRA method can achieve results with higher accuracy, which is compromised by additional effort for adding and removing basis functions to the decomposition.However, DLRA does not require a computationally expensive and memory-intensive offline phase.Furthermore, since DLRA utilizes timedependent basis functions, it exhibits two advantages over POD-Galerkin.First, since basis information can be added and removed in time, the approximation space at a given rank is richer than for POD-Galerkin.This is for example seen when comparing both methods for a fixed rank r = 4 in Fig. 8. DLRA achieves a better approximation quality with fewer degrees of freedom.Second, when the dynamics in the online phase is not captured by the POD-ansatz space, POD-Galerkin requires an expensive retraining while DLRA adapts automatically to such situations.A main advantage of POD-Galerkin, besides the reduced runtime is, that it does not require the derivation of new evolution equations.Moreover, in the case of certain non-linearities, efficient evolution equations for DLRA might not be available.

Square root profile
In this section, we present a realistic test case that implements a water column with an initial square root velocity profile u(0, x, ζ) = √ ζ in the viscous case.Note, that in comparison to Section 7.1 we use an increased slip length and decreased friction coefficient that highlights a stronger coupling of the higher moments to (h, hu m ).The range of parameters is similar to the test cases used in [34], where the relaxation of the model towards different equilibria is investigated.The precise test setup is detailed in Table 3 [42] Table 3: Simulation setup for square root profile test case, compare [42].
As before we sample the solution at two different slip lengths ν = 1 and ν = 100 to set up the POD basis in an offline phase.For each parameter, we sample 800 snapshots, resulting in a total of 1600 snapshots to set up the POD basis as detailed in Section 5.1.
The simulation results are shown in Fig. 13a and Fig. 13b.We see a very good match of the DLRA method with the full HSWME model, while POD-Galerkin and the rHWSME model (with less coefficients r = 4) shows a small deviation from the full model solution.This square root profile example shows that DLRA and its adaptive version are superior for a small number of degrees of freedom.This is especially relevant for realistic scenarios in higher space dimension because here one does not have the luxury of choosing a high number of degrees of freedom due to the memory bottleneck.
Figure 14 also visualizes the difference in the velocity profiles.The same observation as for the macroscopic variables h, u m can be made: DLRA succeeds at reproducing the full model solution of HSWME, while POD-Galerkin and the smaller rHWSME model lead to an error in the velocity profiles.This is true both for a point in the center of the simulation domain in Fig. 14a as well as for a point close to the shock wave in Fig. 14b, where the square root velocity profile is still clearly visible.
Figure 15 shows the speedup obtained depending on the resulting error simulated with the respective method.Each data point represents one choice for the number of the reduced model coefficients r or 10 −10 ≤ ϑ ≤ 10 0 .Note that offline computation costs are neglected here.As expected and shown for the previous test cases, POD-Galerkin has a very fast online computation phase leading to large speedups.The rank-adaptive DLRA can obtain very good accuracy as well even without any offline computation, but at the sacrifice of a more costly online computation.Simply reducing the number of moments in the full model (called rHSWME above) leads to significant errors and is therefore not recommended.
The square root test demonstrates for a more physically relevant velocity profile the power of the rank-adaptive model order reduction via adaptive DLRA.

Summary of numerical results
To summarize and compare the gained insights from numerical test cases, we note that, as expected, reduced moment models (rHSWME) perform well in settings with high friction and smooth initial velocity profiles.This is especially seen in Section 7.1, where rHSWME yields the best performance, achieving a speedup of over 50 compared to the full model to achieve a relative error smaller than 10 −5 .While POD Galerkin demonstrates a speedup of around 20 at the same error level, the additional computational costs required by DLRA lead to a significantly lower speedup of around 5. In situations where the initial velocity profile is irregular or friction is small, POD Galerkin and DLRA exhibit a clear advantage over rHSWME.This is seen in Sections 7.2 and 7.3, where an error level of 10 −5 is surpassed at a speedup of around 50 for POD Galerkin and 16 for DLRA in Section 7.2, and 50 for POD Galerkin and 8 for DLRA in Section 7.3.Overall, POD Galerkin achieves a higher speedup at given error levels compared to DLRA when factoring out offline costs.This is expected, as in this case, POD Galerkin uses a fixed basis, whereas DLRA additionally requires solving basis update equations.However, POD-Galerkin requires additional hyperparameters, while DLRA allows errors to be controlled during computation when using rank-adaptive integrators with a single hyperparameter.Comparing fixed-rank and rank-adaptive DLRA methods, it is observed that an adaptive rank will significantly improve the method's performance.The two factors that benefit the rank-adaptive method are the exact initial conditions of the S-steps and the ability to pick varying ranks over time.

Conclusion
In this work we proposed mass conservative model order reduction methods for the hyperbolic shallow water moment equations that yield fast and accurate solutions.Mass conservation is achieved by decomposing the macroscopic water height and momentum equations from the microscopic higher-order moments and applying model order reduction   For the comparison, we vary the degrees of freedom in the system by adjusting the rank (POD-Galerkin, DLRA) or the number of moments (reduced moments) to r = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 15, 20, 25, 30.For the rank adaptive DLRA, we select 14 different tolerance values on a logarithmic scale from 10 −10 to 10 0 .The relative error denotes the L 2 error of (h, hu m ) at the final time t = 0.05.solely to the microscopic higher-order moment equations.Additionally, the decomposition allows to recover the naive SWE for vanishing truncation rank of the reduced HSWME, resulting in a very accurate reduced model even for a small number of modes.We use two model order reduction methods, namely POD-Galerkin and dynamical low-rank approximations to accelerate the computation and reduce memory footprint.The methods can produce speedups of up to 100 compared to the full HSWME, while the introduced approximation errors are negligible.However, it must be noted, that the HSWME with a reduced number of moments already shows a good performance, if the higher moments vanish quickly.Concluding that DLRA and POD-Galerkin are most efficient if the shearing in the vertical direction is strong and therefore a slow decay in the additional moments is expected.Employing rank-adaptive methods led to further runtime improvements while achieving high accuracy with a smaller number of variables.
The work in this paper opens up possibilities for future work on model reduction for shallow water moment models.An interesting extension would be to use DLRA to generate basis functions for POD.
×r are orthogonal matrices containing the left and right singular vectors, respectively.The latter are also termed modes in the following.According to the Eckart-Young-Mirsky theorem [17, 57] V POD r is the best rank r approximation and the resulting error in the Frobenius norm is rigorously computed from the trailing singular values

Figure 2 :
Figure 2: Sampled snapshot data of the POD, shown for the water height component.The numerical results shown in Fig. 4 compare the numerical results of (1) the full-order HSWME model, (2) the lowest-order SWE model, (3) the new macro-micro decomposition POD-Galerkin, (4) the new macro-micro decomposition DLRA and (5) the HSWME with a reduced number of moments.The full-order HSWME model clearly depicts two waves symmetrically moving left and right, respectively.The simple SWE model uses a constant velocity profile with zero micro structure and thus cannot capture the complex dynamics induced by the bottom friction, which slows down the velocity profile at the bottom.The new macro-micro decomposition POD-Galerkin and DLRA methods, yield optically indistinguishable solutions from the full-order HSWME model, even though only r = 5

Figure
Figure3: Relative approximation error i ∥v(x, t i ) − ṽ(x, t i )∥ L 2 of the snapshot data for space-moment ṽ(x, t) = k α k (t)w k (x) and moment only basis ṽ(x, t) = r k=1 αk (x, t)w k .The online error is given as ∥u(x, t end ) − ũ(x, t end )∥ L 2 , where u = (h, hu m )

Figure 4 :
Figure 4: Water column test case comparison of macroscopic quantities water height h (a) and momentum hu m (b) for full-order HSWME, SWE, POD-Galerkin, and DLRA using rank r = 5.Both reduced models DLRA and POD-Galerkin yield indistinguishable results from the full-order model while the simple SWE model shows insufficient accuracy.

Figure 5 :
Figure 5: Water column test case comparison of velocity profiles close to the domain center (a) and close to the shock wave (b) for full-order HSWME, SWE, POD-Galerkin, and DLRA using rank r = 5.Both reduced models DLRA and POD-Galerkin also yield indistinguishable results from the full-order model while the simple SWE model shows insufficient accuracy and wrong propagation speeds.

Figure 8 :
Figure 8: Smooth wave test case comparison of macroscopic quantities water height h (a) and momentum hu m (b) for full and reduced-moments HSWME, SWE, POD-Galerkin, and DLRA using rank r = 4. Again, the reduced models DLRA and POD-Galerkin yield very good approximations of the full-order model while the simple SWE model fails.

Figure 9 :Figure 10 : 1 − 1 h 1 − 1 1 − 1
Figure 9: Smooth wave test case comparison of velocity profiles close to the domain center (a) and close to the shock wave (b) for full-order HSWME, POD-Galerkin, and DLRA using rank r = 5 (SWE omitted for conciseness).Both reduced models DLRA and POD-Galerkin yield good results in comparison with the full-order model.

Figure 11 : 1 − 1 h
Figure 11: Smooth wave test case time evolution of the total mass ( 1 −1 h dx) relative to its respective initial value plotted for the novel macro-micro decomposition conservative DLRA and a naive (non-conservative) DLRA.The mass is only conserved for the macro-micro decomposition conservative DLRA.

10 − 7 Figure 12 :
Figure 12: Smooth wave test case runtime comparison speedup and error comparison between rHSWME, SWE, DLRA, POD-Galerkin (without offline precomputation phase) in comparison to full HSWME model (100 moments).For the comparison, we vary the degrees of freedom in the system by adjusting the rank (POD-Galerkin, DLRA) or the number of moments (reduced moments) to r = 1, 2, 3, 5, 7, 9, 15, 20, 25, 30.For the rank adaptive DLRA, we select 14 different tolerance values on a logarithmic scale from 10 −10 to 10 0 .The relative error denotes the L 2 error of u = (h, hu m ) at the final time t = 0.2.

Figure 13 :
Figure 13: Square root test case comparison of macroscopic variables h (a) and hu m (b) for full-order HSWME, SWE, POD-Galerkin, and DLRA using rank r = 4 and a low-order rHSWME model with only 4 moments.DLRA shows excellent accuracy and POD-Galerkin yield good results in comparison with the full-order model.

Figure 14 :
Figure 14: Square root test case comparison of velocity profiles close to the domain center (a) and close to the shock wave (b) for full-order HSWME, POD-Galerkin, DLRA using r = 4 and a low-order rHSWME model with only r = 4 modes (SWE omitted for conciseness).Again DLRA shows excellent accuracy and POD-Galerkin yield good results in comparison with the full-order model.
Julian Koellermeier: initial idea, methodology, implementation of HSWME, setup of test cases Philipp Krah: implementation of POD-Galerkin, simulation/setup of numerical tests, visualization Jonas Kusch: initial idea, methodology, implementation of HSWME, implementation of DLRA, simulation/setup of numerical tests, visualization