# Momentum-conserving ROMs for the incompressible Navier-Stokes equations

Henrik K. E. Rosenberger1 and Benjamin Sanderse2 1 Centrum Wiskunde &
Informatica
Science Park 123, Amsterdam, The Netherlands
email<EMAIL_ADDRESS>www.cwi.nl/people/henrik-rosenberger and 2
Centrum Wiskunde & Informatica
Science Park 123, Amsterdam, The Netherlands
email<EMAIL_ADDRESS>www.thinkingslow.nl

###### Abstract

Projection-based model order reduction of an ordinary differential equation
(ODE) results in a projected ODE. Based on this ODE, an existing reduced-order
model (ROM) for finite volume discretizations satisfies the underlying
conservation law over arbitrarily chosen subdomains. However, this ROM does
not satisfy the projected ODE exactly but introduces an additional
perturbation term. In this work, we propose a novel ROM with the same
subdomain conservation properties which indeed satisfies the projected ODE
exactly.

We apply this ROM to the incompressible Navier-Stokes equations and show with
regard to the mass equation how the novel ROM can be constructed to satisfy
algebraic constraints.

Furthermore, we show that the resulting mass-conserving ROM allows us to
derive kinetic energy conservation and consequently nonlinear stability, which
was not possible for the existing ROM due to the presence of the perturbation
term.

###### keywords:

Model order reduction, finite volume method, structure preservation,
incompressible Navier-Stokes equations

Henrik K. E. Rosenberger and Benjamin Sanderse

##  1 INTRODUCTION

For many applications such as optimization, uncertainty quantification and
real-time simulation of large-scale systems, classical numerical methods such
as finite elements or finite volumes are prohibitively expensive. To mitigate
the computational costs, several techniques have been developed to reduce the
complexity of the full order models (FOM) developed via such classical
methods, yielding so-called reduced order models (ROM). Many of these
techniques are projection-based, i.e., the FOM is projected onto a lower-
dimensional subspace spanned by a reduced basis. Such a reduced basis can be
obtained, e.g., via Proper Orthogonal Decomposition of FOM snapshots [1].

While ROMs have been developed successfully for elliptic and parabolic systems
of partial differential equations, deriving stable ROMs for hyperbolic systems
is challenging [2]. Stability problems of these ROMs are often ascribed to
lacking structure preservation. While many FOMs feature physical structures of
the underlying systems, many model reduction techniques do not preserve such
structures [2].

Therefore, many approaches to obtain stable reduced order models refine model
reduction techniques such that they do preserve such physical structures [3].
Peng et al. [4] propose a POD-Galerkin ROM for Hamiltonian systems that
preserve the systems’ symplecticity. Chan [5] constructs ROMs for nonlinear
conservation laws that exhibit a semi-discrete entropy dissipation. Sanderse
[6] achieves nonlinear stability of a POD-Galerkin ROM for the incompressible
Navier-Stokes equations by preserving the kinetic energy evolution and
additionally conserves the global momentum.

Carlberg et al. [7] address another aspect of structure preservation specific
to finite volume discretizations. By construction, these discretizations
satisfy the underlying conservation law over all finite volumes. Projection-
based ROMs, on the other hand, satisfy this conservation law generally over
none of the finite volumes. While it is in most cases impossible to define
ROMs that satisfy the conservation law over all finite volumes, Carlberg et
al. propose ROMs that satisfy the conservation law at least over a few
subdomains. These ROMs are formulated via a constrained optimization problem
(COP) and are shown to be more accurate than ROMs without subdomain
conservation. Furthermore, the COP framework has been extended by Schein et
al. [8] to include other types of constraints.

In this work, we combine the concept of conservation over subdomain in [7] and
the nonlinear stable model reduction in [6]. Unfortunately, the subdomain
conservative ROM in [7] is not suitable for this purpose. As shown in [7], the
solution to the COP can equivalently be expressed by an ODE. This ODE,
however, differs from the conventional Galerkin ROM ODE by a perturbation
term. This perturbation term impedes the derivation of nonlinear stability as
described in [6]. Therefore, we propose a novel subdomain conservative ROM
which does not introduce any perturbation terms.

We apply the novel ROM to the energy-conserving finite volume discretization
of the incompressible Navier-Stokes equations proposed in [6]. Regarding the
momentum equation, we can directly apply the novel ROM. Regarding the mass
equation, however, we have to extend the novel ROM to treat algebraic
constraints. The resulting mass-conserving ROM is shown to mimick the kinetic
energy evolution of the FOM. This evolution states that the global kinetic
energy does not increase over time, hence implying nonlinear stability of the
ROM.

This article is structured as follows. First, we summarize the subdomain
conservative ROM proposed in [7] in Section 2. Then, we derive our novel
subdomain conservative ROM in Section 3, in the context of a general
conservation law. In Section 4, we apply our ROM to the energy-conserving
finite volume discretization of the incompressible Navier-Stokes equations
proposed in [6] and show that the ROM preserves the kinetic energy evolution
of the FOM. We summarize our findings in Section 5 and give an outlook on
possible uses of the proposed ROMs in Section 6.

##  2 SUBDOMAIN CONSERVATION FOR GENERAL CONSERVATION LAWS

In this section, we summarize the concept of subdomain conservation and the
subdomain conservative ROM described in [7]. For convenience only, we use a
simplified notation and consider only scalar-valued conservation laws.

###  2.1 Subdomain conservation for finite volume discretizations

In the understanding of Carlberg et al. [7], a scalar-valued quantity $u$ is
conserved over a domain $\omega$ if it satisfies the integral form of a
conservation law over this domain

$\displaystyle\frac{\textrm{d}}{\textrm{d}t}\int_{\omega}u\textrm{\;d}V+\int_{\partial\omega}g\cdot
n\textrm{\;d}S=\int_{\omega}s\textrm{\;d}V.$ (1)

where $g=g(u)$ describes the flux across the boundary and $s$ is a term
describing sources and sinks.

For example, a quantity $u$ satisfying the differential form of a conservation
law

$\displaystyle\frac{\partial}{\partial t}u+\nabla\cdot g=s\hskip
20.00003pt\textrm{ in }\Omega,$ (2)

is conservative over any subdomain of $\Omega$, because (2) implies (1) for
all $\omega\subset\Omega$.

Let us now consider a finite volume discretization of (2) on a mesh
$\mathcal{M}=\\{\Omega_{j}\\}_{j=1}^{N_{\mathcal{M}}}$ that completely covers
$\Omega$. Then, we can define the time-dependent state-vector
$\mathbf{u_{h}}=\begin{bmatrix}u_{1}\;\dots\;u_{j}\;\dots\;u_{N_{\mathcal{M}}}\end{bmatrix}^{T}$
of cell averages

$\displaystyle
u_{j}=\frac{1}{\left|\Omega_{j}\right|}\int_{\Omega_{j}}u\textrm{\;d}V\hskip
10.00002ptj=1,\dots,N_{\mathcal{M}},$ (3)

and the time-dependent vector
$\mathbf{f}=\begin{bmatrix}-f^{g}_{1}+f^{s}_{1},\;\dots\;-f^{g}_{j}+f^{s}_{j},\;\dots\;-f^{g}_{N_{\mathcal{M}}}+f^{s}_{N_{\mathcal{M}}}\end{bmatrix}^{T}$
consisting of the net flux

$\displaystyle f^{g}_{j}$
$\displaystyle=\frac{1}{\left|\Omega_{j}\right|}\int_{\partial\Omega_{j}}g^{FV}\cdot
n\textrm{\;d}S,$ (4) and the net sinks and sources $\displaystyle f^{s}_{j}$
$\displaystyle=\frac{1}{\left|\Omega_{j}\right|}\int_{\Omega_{j}}s^{FV}\textrm{\;d}V,$
(5)

of each cell, divided by the respective cell-volume. Here, $g^{FV}$ and
$s^{FV}$ are approximations of $g$ and $s$ in (2).

Together, the state vector $\mathbf{u_{h}}(t)$ and the RHS vector
$\mathbf{f}(\mathbf{u_{h}},t)$ form the ordinary differential equation (ODE)
system

$\displaystyle\frac{\textrm{d}}{\textrm{d}t}\mathbf{u_{h}}=\mathbf{f}(\mathbf{u_{h}},t),$
(6)

which we can interpret as a system of conservation laws in integral form over
all finite volume cells.

Even more, the vector $\mathbf{u_{h}}$ represents a function on $\Omega$ which
is constant over each finite volume. This function is subdomain conservative
over all finite volume cells with respect to the approximated fluxes $g^{FV}$,
and sink and source terms $s^{FV}$. Therefore, we also denote $\mathbf{u_{h}}$
itself as subdomain conservative over all finite volume cells.

In the following, it will be useful to express the ODE (6) equivalently as

$\displaystyle\mathbf{r}\left(\frac{\textrm{d}}{\textrm{d}t}\mathbf{u_{h}},\mathbf{u_{h}},t\right)=0,$
(7)

with the residual
$\mathbf{r}(\mathbf{v},\mathbf{w},t)=\mathbf{v}-\mathbf{f}(\mathbf{w},t)$.

###  2.2 Subdomain conservation for Galerkin projection ROMs

Let us now investigate the subdomain conservation properties of a model ROM
approximating the FOM (6). Given an orthonormal basis
$\Phi\in\mathbb{R}^{N_{\mathcal{M}}\times p},\;p\ll N_{\mathcal{M}}$, we
consider the ROM

$\displaystyle\mathbf{u_{r}}(t)=\Phi\mathbf{a}(t)\approx\hskip
10.00002pt\mathbf{u_{h}}(t),$ (8)

that satisfies the Galerkin projection of the FOM ODE (6)

$\displaystyle\Phi^{T}\frac{\textrm{\;d}}{\textrm{\;d}t}\mathbf{u_{r}}=\frac{\textrm{\;d}}{\textrm{\;d}t}\mathbf{a}=\Phi^{T}\mathbf{f}(\Phi\mathbf{a},t)\hskip
10.00002pt\Leftrightarrow\hskip
10.00002pt\Phi^{T}\mathbf{r}\left(\Phi\frac{\textrm{d}}{\textrm{d}t}\mathbf{a},\Phi\mathbf{a},t\right)=0.$
(9)

As shown in [9], this ODE system, equipped with some initial condition
$\mathbf{a}(0)=\mathbf{a}_{0}$, can be equivalently described by the
optimization problem

$\displaystyle\frac{\textrm{d}}{\textrm{d}t}\mathbf{a}=\arg\min_{\mathbf{b}\in\mathbb{R}^{p}}\|r\left(\Phi\mathbf{b},\Phi\mathbf{a},t\right)\|_{2},$
(10)

and the same initial condition. Hence, the ROM $\mathbf{u_{r}}(t)$ minimizes
the violation of the cell-wise conservation laws (6) but this minimum is not
guaranteed to vanish. Consequently, the ROM $\mathbf{u_{r}}(t)$ is not
guaranteed to be subdomain conservative over any finite volume.

Hence, instead of aiming at subdomain conservation over all finite volumes,
Carlberg et al. propose to require conservation over a set of a few subdomains
each consisting of one or more finite volumes.

To this end, we consider a decomposed mesh $\mathrm{bar}\mathcal{M}$ of
$N_{\mathrm{bar}\mathcal{M}}$ arbitrary subdomains
$\mathrm{bar}\Omega_{k}=\cup_{j\in S_{k}}\Omega_{j}$ of finite volumes
$\Omega_{j}$ where $S_{k}\subset\\{1,\dots,N_{\mathcal{M}}\\}$. These
subdomains $\mathrm{bar}\Omega_{h}$ can overlap and do not need to be
connected. An example of such a decomposed mesh is depicted in Fig. 1.

Figure 1: Example of three subdomains.

On this grid, we define, similarly to the finite volume quantities (3), (4)
and (5), a time-dependent state vector
$\mathbf{\mathrm{bar}u_{h}}=\begin{bmatrix}\mathrm{bar}u_{1}\;\dots\;\mathrm{bar}u_{k}\;\dots\;\mathrm{bar}u_{N_{\mathrm{bar}\mathcal{M}}}\end{bmatrix}^{T}$
of subdomain averages

$\displaystyle\mathrm{bar}u_{k}=\frac{1}{\left|\mathrm{bar}\Omega_{k}\right|}\int_{\mathrm{bar}\Omega_{k}}u\textrm{\;d}V=\frac{1}{\left|\mathrm{bar}\Omega_{k}\right|}\sum_{j\in
S_{k}}\int_{\Omega_{j}}u\textrm{\;d}V=\frac{1}{\left|\mathrm{bar}\Omega_{k}\right|}\sum_{j\in
S_{k}}\left|\Omega_{j}\right|u_{j},$ (11)

and analogously a time-dependent RHS vector $\mathbf{\mathrm{bar}f},$
consisting of the net flux and the net sinks and sources of each subdomain,
divided by the respective subdomain volume. We can write
$\mathbf{\mathrm{bar}u_{h}}$ and $\mathbf{\mathrm{bar}f}$ in matrix vector
notation as

$\displaystyle\mathbf{\mathrm{bar}u_{h}}$ $\displaystyle=C^{T}\mathbf{u_{h}}$
$\displaystyle\mathbf{\mathrm{bar}f}$ $\displaystyle=C^{T}\mathbf{f},$ (12)

for some $C\in\mathbb{R}^{N_{\mathcal{M}}\times N_{\mathrm{bar}\mathcal{M}}}$.

As a result, subdomain conservation over all subdomains of
$\mathrm{bar}\mathcal{M}$ is described by the $C^{T}$-premultiplied ODE system

$\displaystyle C^{T}\frac{d}{dt}\mathbf{u_{h}}=C^{T}\,\mathbf{f}\hskip
10.00002pt\Leftrightarrow\hskip
10.00002ptC^{T}\,\mathbf{r}\left(\frac{d}{dt}\mathbf{u_{h}},\mathbf{u_{h}},t\right)=0.$
(13)

Based on the Galerkin projection ROM (9), we want to construct a ROM that also
satisfies the subdomain conservation constraint (13). In the following
section, we summarize the subdomain conservative ROM proposed in [7]. Our
novel subdomain conservative ROM is presented in Section 3.

###  2.3 Existing approach for subdomain conservative ROMs: constrained
optimization problem formulation

Carlberg et al. [7] suggest to use the interpretation of the ROM ODE (9) as
minimization problem (10). Adding the conservation constraint (13) to this
optimization problem, we get the constrained optimization problem

$\displaystyle\text{minimize}_{\,\mathbf{b}\in\mathbb{R}^{p}}$
$\displaystyle\|r(\Phi\mathbf{b},\Phi\mathbf{a},t)$ $\displaystyle\|_{2}$ (14)
subject to $\displaystyle C^{T}\mathbf{r}(\Phi\mathbf{b},\Phi\mathbf{a},t)$
$\displaystyle=0.$ (15)

A sufficient condition for feasibility of this COP is given in [7, Proposition
5.1]. If the COP is feasible, the solution is equivalently expressed by the
ODE

$\displaystyle\frac{\textrm{d}}{\textrm{d}t}\mathbf{a}=\Phi^{T}\mathbf{f}(\Phi\mathbf{a},t)+\mathbf{f^{*}}(\Phi\mathbf{a},t)$
(16)

with the perturbation term

$\displaystyle\mathbf{f}^{*}(\Phi\mathbf{a},t)=(C^{T}\Phi)^{+}[C^{T}-C^{T}\Phi\Phi^{T}]\mathbf{f}(\Phi\mathbf{a},t),$
(17)

and the Moore-Penrose inverse $(C^{T}\Phi)^{+}$ of $C^{T}\Phi$.

Except for the perturbation term, this ODE is equivalent to the conventional
Galerkin ROM ODE (9). The perturbation term, however, lacks a clear physical
motivation and impedes the preservation of the kinetic energy evolution as we
will see in Section 4.4. Therefore, we propose a subdomain conservative ROM
without any perturbation terms in the next section.

##  3 NOVEL APPROACH FOR SUBDOMAIN CONSERVATIVE ROMS: BASIS MODIFICATION

###  3.1 Novel approach

The basic idea of our novel approach is to merge the ROM ODE (9) and the
conservation constraint (13) in one linear system for
$\frac{\textrm{d}}{\textrm{d}t}\mathbf{u_{r}}$,

$\displaystyle\left[\Phi\;\;C\right]^{T}\frac{\textrm{d}}{\textrm{d}t}\mathbf{u_{r}}=\left[\Phi\;\;C\right]^{T}\mathbf{f}(\mathbf{u_{r}},t).$
(18)

To satisfy this linear system, we define an orthogonal basis
$\tilde{\Phi}\,\in\mathbb{R}^{N_{\mathcal{M}}\times
H},\;H:=\operatorname{rank}([\Phi\;\;C])$ that spans the same linear subspace
as $[\Phi\;\;C]$ and equip our novel ROM with this basis,

$\displaystyle\mathbf{u_{r}}(t)=\tilde{\Phi}\mathbf{\tilde{a}}(t),\;\mathbf{\tilde{a}}(t)\in\mathbb{R}^{H}.$
(19)

Then, the Galerkin ROM ODE

$\displaystyle\tilde{\Phi}^{T}\frac{\textrm{d}}{\textrm{d}t}\tilde{\Phi}\mathbf{\tilde{a}}=\frac{\textrm{d}}{\textrm{d}t}\mathbf{\tilde{a}}=\tilde{\Phi}^{T}f(\tilde{\Phi}\mathbf{\tilde{a}},t),$
(20)

satisfies the system (18). Hence, this ROM is subdomain conservative without
introducing any perturbation terms.

###  3.2 View on POD bases

A popular method to construct ROM bases is POD, which exhibits an optimality
property: given a snapshot matrix
$X=[\mathbf{u_{h}^{0}}\;\;\dots\;\;\mathbf{u_{h}^{K-1}}]\in\mathbb{R}^{N_{\mathcal{M}}\times
K}$, a POD basis $\Phi_{\mathrm{POD}}\in\mathbb{R}^{N_{\mathcal{M}}\times p}$
is known to minimize the sum of best approximation errors of the snapshots in
$X$ among all orthogonal matrices in $\mathbb{R}^{N_{\mathcal{M}}\times p}$
[1].

On the other hand, the key element of our novel approach is the modification
of the ROM basis in consideration of the predefined subdomains. In order that
the Galerkin ROM ODE (20) includes the subdomain conservation constraint (13),
the span of the basis $\tilde{\Phi}$ must include the span of $C$. This
requirement can be expressed as the equivalence of $C$ and the orthogonal
projection of $C$ onto $\tilde{\Phi}$, $C=\tilde{\Phi}\tilde{\Phi}^{T}C$.
Combining the minimization idea of POD and the subdomain conservation
constraint, we find the constrained optimization problem for the ROM basis
$\tilde{\Phi}$

$\displaystyle\text{minimize}_{\;\Xi\in\mathbb{R}^{N_{\mathcal{M}}\times q}}$
$\displaystyle\sum_{j=0}^{K-1}\|\mathbf{u_{h}^{j}}-\Xi\Xi^{T}\mathbf{u_{h}^{j}}$
$\displaystyle\|^{2}$ (21) subject to $\displaystyle\Xi^{T}\Xi$
$\displaystyle=I$ (22) and $\displaystyle C$ $\displaystyle=\Xi\Xi^{T}C.$ (23)

The orthogonal basis that spans the same subspace as
$[\Phi_{\mathrm{POD}}\;\;C]$ does generally not solve this constrained
optimization problem.

A solution to this constrained optimization problem is proposed in [10]. This
solution is defined as

$\displaystyle\tilde{\Phi}=[Q_{1}\;\;Q_{2}V],$ (24)

where $Q_{1}\in\mathbb{R}^{N_{\mathcal{M}}\times H_{C}}$ and
$Q_{2}\in\mathbb{R}^{N_{\mathcal{M}}\times(N_{\mathcal{M}}-H_{C}})$ with
$H_{C}=\operatorname{rank}(C)$ are obtained from the QR decomposition of
$C^{T}$,

$\displaystyle C=[Q_{1}\;\;Q_{2}]\begin{bmatrix}R_{1}\\\
0\end{bmatrix},\;R_{1}\in\mathbb{R}^{H_{C}\times N_{\mathrm{bar}\mathcal{M}}}$
(25)

and $V\in\mathbb{R}^{(N_{\mathcal{M}}-H_{C})\times(q-H_{C})}$ consists of the
first $q-H_{C}$ modes of the POD basis of $Q_{2}^{T}X$.

In view of Section 3.1, $V$ can be interpreted as the initial basis $\Phi$.
Observe that in contrast to $\Phi_{\mathrm{POD}}$, the matrix $V$ depends on
$C$.

###  3.3 View on least-squeares Petrov-Galerkin ROMs

In [7], Carlberg et al. discuss besides the Galerkin ROM also a least-squares
Petrov-Galerkin (LSPG) ROM. The two ROMs differ in the order of performing the
ROM approximation and the time discretization in their derivations. While the
Galerkin ROM is obtained by introducing the ROM approximation on the time-
continuous model and discretizing in time afterwards, the LSPG ROM is obtained
by first discretizing in time and then performing the ROM approximation on the
time-discrete model.

So far, we have introduced our novel approach only for the Galerkin ROM. In
fact, our approach does not directly work in the context of the LSPG ROM. For
the LSPG ROM trial and test basis are not the same. As a consequence, the
choice of the basis $\eqref{eq:def tilde phi}$ as trial basis does generally
not imply that the span of the test basis includes the span of $C$. Hence, the
subdomain conservation constraint (13) is generally not satisfied.

###  3.4 Subdomain conservation for invariants

We want to highlight a property of our novel ROM for the special case of
invariants. Invariants are quantities that are constant over time, e.g., the
global integral of a conservation quantity with periodic boundary conditions
and in absence of external forces. If constraints in the constraint matrix $C$
describe such invariants, then the coefficients in $\mathbf{\tilde{a}}(t)$ in
the expansion (19) that correspond to those invariant constraints are
constant. Hence, we do not need to integrate these coefficients via the ODE
(20), but can compute their values from the initial conditions. This
observation can be seen as a physical justification for a ROM Ansatz with a
constant offset,
$\mathbf{u_{r}}(t)=\tilde{\Phi}\mathbf{\tilde{a}}(t)+\mathbf{u_{r}^{0}}$.

##  4 MOMENTUM CONSERVATION OVER SUBDOMAINS FOR THE INCOMPRESSIBLE NAVIER-
STOKES EQUATIONS

###  4.1 Introduction to the FOM

We consider the energy-conserving finite volume discretization of the
incompressible Navier-Stokes equations with periodic boundary conditions
proposed in [6],

$\displaystyle M_{h}\mathbf{V_{h}}(t)$ $\displaystyle=0,$ (26)
$\displaystyle\Omega_{h}\frac{\textrm{d}}{\textrm{d}t}\mathbf{V_{h}}(t)$
$\displaystyle=\mathbf{F_{h}^{CD}}(\mathbf{V_{h}}(t),t)-G_{h}\mathbf{p_{h}}(t).$
(27)

The vectors $\mathbf{V_{h}}(t)\in\mathbb{R}^{N_{V}}$ and
$\mathbf{p_{h}}\in\mathbb{R}^{N_{p}}$ describe the velocity and the pressure,
respectively. The matrices $M_{h}\in\mathbb{R}^{N_{p}\times N_{V}}$ and
$G_{h}\in\mathbb{R}^{N_{V}\times N_{p}}$ are the discretizations of the
divergence and the gradient operator, respectively, and satisfy the duality
property $M_{h}=-G_{h}^{T}$. The matrix $\Omega_{h}\in\mathbb{R}^{N_{V}\times
N_{V}}$ is diagonal with the sizes of the finite volumes on its diagonal, so
symmetric positive definite. The term
$\mathbf{F_{h}^{CD}}(\mathbf{V_{h}}(t),t)\in\mathbb{R}^{N_{V}}$ comprises
convective, diffusive and body force contributions.

While the momentum equation (27) is an ODE, the mass equation (26) is an
algebraic equation. The subdomain conservation ROM method described in Section
3, however, addresses conservation laws of ODE structure only. Therefore, we
apply the subdomain conservation method only to the momentum equation and deal
with the mass equation separately in Section 4.3.

###  4.2 Momentum conservation over subdomains

To construct a ROM that conserves momentum over a set of predefined subdomains
represented by the matrix $C\in\mathbb{R}^{N_{V}\times
N_{\hat{\mathcal{M}}}}$, we define the approximation of the velocity

$\displaystyle\mathbf{V_{r}}(t)=\phi\mathbf{a}(t)\hskip
10.00002pt\approx\mathbf{V_{r}}(t),\hskip
10.00002pt\phi\in\mathbb{R}^{N_{V}\times R_{V}},\hskip
10.00002pt\mathbf{a}(t)\in\mathbb{R}^{R_{V}},\hskip 10.00002ptR_{V}\ll N_{V}.$
(28)

In contrast to Section 3, where we introduced the orthogonal basis
$\tilde{\Phi}$ with respect to the $L^{2}$-inner product, we construct the
velocity basis $\phi$ to be orthogonal with respect to the
$\Omega_{h}$-weighted inner product, i.e. $\phi^{T}\Omega_{h}\phi=I$. As a
result, the definition of $\phi$,

$\displaystyle\phi=[\tilde{Q}_{1}\;\;W],$ (29)

differs slightly from the definition of $\tilde{\Phi}$ (24). Here,
$\tilde{Q}_{1}=\Omega_{h}^{-1/2}Q_{1}$ with the QR decomposition of
$\Omega_{h}^{-1/2}C$,

$\displaystyle\Omega_{h}^{-1/2}C=[Q_{1}\;\;Q_{2}]\begin{bmatrix}R_{1}\\\
0\end{bmatrix}=Q_{1}R_{1},$ (30)

and $W$ consists of the first $R_{V}-\operatorname{rank}(C)$ columns of the
$\Omega_{h}$-orthogonal POD basis $\tilde{U}$ of a snapshot matrix
$X=[\mathbf{V_{h}^{0}}\;\;\dots\;\;\mathbf{V_{h}^{K-1}}]\in\mathbb{R}^{N_{V}\times
K}$ that is $\Omega_{h}$-orthogonal to $\tilde{Q}_{1}$. In detail:

We project the snapshot matrix $X$ onto the subspace $\Omega_{h}$-orthogonal
to $\tilde{Q}_{1}$: $\displaystyle\tilde{X}$
$\displaystyle=[I-\tilde{Q}_{1}\tilde{Q}_{1}^{T}\Omega_{h}]X.$ (31) We
transform to include the $\Omega_{h}$-inner product: $\displaystyle\hat{X}$
$\displaystyle=\Omega_{h}^{1/2}\tilde{X}.$ (32) We compute the SVD of
$\hat{X}$: $\displaystyle\hat{X}$
$\displaystyle=\hat{U}\hat{\Sigma}\hat{V}^{T}.$ (33) We transform back:
$\displaystyle\tilde{U}$ $\displaystyle=\Omega_{h}^{-1/2}\hat{U}.$ (34)
Finally, we truncate: $\displaystyle W$
$\displaystyle=[\tilde{U}]_{R_{V}-\operatorname{rank}(C^{T})}.$ (35)

This construction is a generalization of the construction in [6, Appendix C].

By inserting the velocity approximation (28) into the momentum equation (27)
and Galerkin-projecting, we find the ODE

$\displaystyle\phi^{T}\Omega_{h}\frac{\textrm{d}}{\textrm{d}t}\phi\mathbf{a}(t)=\frac{\textrm{d}}{\textrm{d}t}\mathbf{a}(t)=\phi^{T}\mathbf{F_{h}^{CD}}(\phi\mathbf{a}(t),t)-\phi^{T}G_{h}\mathbf{p}_{h}(t).$
(36)

The right-hand side of this ODE contains a pressure term. We will eliminate
this pressure term in the next section as a by-product of enforcing mass
conservation.

###  4.3 Mass conservation

Inserting the velocity approximation (28) into the mass equation (26), we find

$\displaystyle M_{h}\phi\mathbf{a}(t)=0.$ (37)

To satisfy mass conservation, we observe that the ROM (28) satisfies the mass
equation (37) for all $\mathbf{a}(t)\in\mathbb{R}^{R_{V}}$, if the ROM basis
$\phi$ is divergence-free, i.e. $M_{h}\phi=0$.

To exploit this insight if $\phi$ is not divergence-free, we decompose the
basis into a divergence-free component $\phi_{0}$ and a basis orthogonal to
the space of divergence-free vectors, $\phi_{\bot}$. For this purpose, we
compute the QR decomposition of $(M_{h}\phi)^{T}$,

$\displaystyle(M_{h}\phi)^{T}=[Q_{1}^{M}\;\;Q_{2}^{M}]\begin{bmatrix}R_{1}^{M}\\\
0\end{bmatrix}=Q_{1}^{M}R_{1}^{M}$ (38)

and decompose

$\displaystyle\mathbf{a}(t)=Q_{1}^{M}\mathbf{a_{1}}(t)+Q_{2}^{M}\mathbf{a_{2}}(t).$
(39)

Then, we find
$\phi\mathbf{a}(t)=\phi_{0}\mathbf{a_{2}}(t)+\phi_{\bot}\mathbf{a_{1}}(t)$
with the divergence-free basis $\phi_{0}=\phi Q_{2}^{M}$ and the basis
$\phi_{\bot}=\phi Q_{1}^{M}$ orthogonal to the divergence-free subspace.
Inserting the decomposition (39) in the mass equation (37), we find

$\displaystyle
M_{h}\phi\mathbf{a}(t)=(R_{1}^{M})^{T}(Q_{1}^{M})^{T}[Q_{1}^{M}\mathbf{a_{1}}(t)+Q_{2}^{M}\mathbf{a_{2}}(t)]=(R_{1}^{M})^{T}\mathbf{a_{1}}(t).$
(40)

Since $(R_{1}^{M})^{T}$ has full column rank, the velocity approximation (28)
satisfies the mass equation (37), if and only if we set $\mathbf{a_{1}}(t)=0$.
As a result, the velocity approximation simplifies to $V_{r}(t)=\phi
Q_{2}^{M}\mathbf{a_{2}}(t)=\phi_{0}\mathbf{a_{2}}(t)$. The remaining
coefficients $\mathbf{a_{2}}(t)$ are computed via the
$(Q_{2}^{M})^{T}$-premultiplied ODE (36),

$\displaystyle(Q_{2}^{M})^{T}\frac{\textrm{d}}{\textrm{d}t}Q_{2}^{M}\mathbf{a_{2}}(t)=\frac{\textrm{d}}{\textrm{d}t}\mathbf{a_{2}}(t)=(Q_{2}^{M})^{T}\phi^{T}\mathbf{F_{h}^{CD}}(\phi_{0}\mathbf{a_{2}}(t),t).$
(41)

The pressure term is omitted because
$(Q_{2}^{M})^{T}\phi^{T}G_{h}=\phi_{0}^{T}G_{h}=(M_{h}\phi_{0})^{T}=0$.

To solve this ODE, we propose to employ an energy-conserving Runge-Kutta
method, e.g., [11].

###  4.4 Kinetic energy conservation

As in [6], we define the kinetic energy of the ROM as

$\displaystyle
K_{r}(t)=\frac{1}{2}\mathbf{V_{r}}(t)^{T}\Omega_{h}\mathbf{V_{r}}(t)=\frac{1}{2}\|\mathbf{a_{2}}(t)\|_{2}^{2}.$
(42)

Using the ODE (41), we calculate

$\displaystyle\frac{\textrm{d}}{\textrm{d}t}K_{r}(t)$
$\displaystyle=\mathbf{a_{2}}(t)^{T}\frac{\textrm{d}}{\textrm{d}t}\mathbf{a_{2}}(t)=\mathbf{a_{2}}(t)^{T}\phi_{0}^{T}\mathbf{F_{h}^{CD}}(\phi_{0}\mathbf{a_{2}}(t),t),$
(43)

where there pressure term is eliminated due to the mass equation (37).

As shown in [6], the right-hand side evaluates in absence of body forces to
$-\nu\|Q_{r}\mathbf{a_{2}}\|_{2}^{2}$ with some matrix
$Q_{r}\in\mathbb{R}^{N_{V}\times R_{V}}$. Hence, the kinetic energy is
constant in the inviscid limit and decreases otherwise. Consequently, the ROM
is nonlinearly stable.

In contrast, for the ROM proposed in [7], the kinetic energy evolution is
given by

$\displaystyle\frac{\textrm{d}}{\textrm{d}t}K_{r}(t)=-\nu\|Q_{r}\mathbf{a_{2}}(t)\|_{2}^{2}+\mathbf{a_{2}}(t)^{T}\mathbf{F_{h}^{CD,*}}(\phi\mathbf{a_{2}}(t),t),$
(44)

provided that we use the approach described in Section 4.3 to enforce mass
conservation. Because of the perturbation term
$\mathbf{F_{h}^{CD,*}}(\phi\mathbf{a_{2}}(t),t)$ corresponding to (17), it is
unclear whether the global kinetic energy increases over time and and whether
nonlinear stability can be inferred.

##  5 CONCLUSION

In this article, we have proposed a novel approach to develop reduced order
models (ROMs) for finite volume discretizations that preserve the underlying
conservation law over arbitrary subdomains. The key element of this approach
is the modification of the ROM basis in consideration of the predefined
subdomains. This modification incorporates the subdomain conservation
constraints into the test basis and thereby enforces their satisfaction. The
main advantage of our proposed ROM is that it does not introduce perturbation
terms and can be simply defined via the conventional Galerkin ROM ODE.

For the special case that conservation quantities are constant in time over
the chosen subdomain, we have shown that the ROMs can be computed more
efficiently by introducing a physically motivated constant offset.

Furthermore, we have applied the novel subdomain conservation approach to an
energy-conserving finite volume discretization of the incompressible Navier-
Stokes equations. To enforce the mass equation on the ROM level, we have
proposed a method based on a QR decomposition (which can also be generalized
to other algebraic constraints). From the mass conservation, we have inferred
the kinetic energy evolution which implies nonlinear stability of the ROM.
This makes it stand out from existing subdomain conservative ROMs that involve
perturbation terms, for which such an argument probably does not hold.

##  6 OUTLOOK

In this paper we have not addressed the choice of subdomains and their effect
on the accuracy of the resulting ROMs. We believe that subdomain conservation
could improve the generalization accuracy of ROMs. In particular, we see a
potential in choosing a small number of subdomains each consisting of a single
finite volume. Consequently, the ROM would be simulated in the classical
finite volume sense over these finite volumes, while the remaining finite
volumes are simulated in the Galerkin ROM fashion. This way, we can decompose
the simulation domain into areas where we have confidence in the reduced basis
and areas where we foresee that a high-fidelity method is required.

This interpretation can be seen as a generalization of heterogeneously refined
grids in classical numerical methods. Such heterogeneous grids have highly
refined grid cells where they are believed to be necessary and coarse grid
cells everywhere else. The size of the grid cells is either computed apriori
based on prior knowledge or adaptively based on a posteriori error estimates.
In a similar fashion we could determine the conservation subdomains based on
prior knowledge, or based on the analysis of snapshot data of the FOM
simulation.

An example of informative prior knowledge is a parametrized external force
that acts locally in the domain. In the area where this force acts, the
solution likely changes significantly over variation of the parameters. Hence,
snapshot data of a given set of parametrizations might not generalize well to
other parametrizations. Therefore, we suggest to compute the ROM via finite
volumes in this area.

As an idea to determine conservation subdomains based on snapshot data, we
suggest to compute in which finite volumes the conservation quantity changes
the most over time. We interpret the conservation laws over these finite
volumes as the most important ones to enforce and hence suggest to choose
these finite volumes as conservation subdomains.

Apart from these ideas to improve the global accuracy, enforcing conservation
over single finite volumes could be used to improve the ROM accuracy locally.
This approach could be useful for applications that are characterized by a
particular interest in simulation results only in a small portion of the
computation domain, e.g., flow around an airfoil.

## References

  * [1] S. Volkwein, Model reduction using proper orthogonal decomposition, Lecture Notes, Institute of Mathematics and Scientific Computing, University of Graz. see http://www. uni-graz. at/imawww/volkwein/POD. pdf 1025 (2011).
  * [2] B. M. Afkham, N. Ripamonti, Q. Wang, J. S. Hesthaven, Conservative model order reduction for fluid flow, in: Quantification of Uncertainty: Improving Efficiency and Technology, Springer, 2020, pp. 67–99.
  * [3] J. S. Hesthaven, C. Pagliantini, G. Rozza, Reduced basis methods for time-dependent problems, Acta Numerica 31 (2022) 265–345.
  * [4] L. Peng, K. Mohseni, Symplectic model reduction of Hamiltonian systems, SIAM Journal on Scientific Computing 38 (1) (2016) A1–A27.
  * [5] J. Chan, Entropy stable reduced order modeling of nonlinear conservation laws, Journal of Computational Physics 423 (2020) 109789.
  * [6] B. Sanderse, Non-linearly stable reduced-order models for incompressible flow with energy-conserving finite volume methods, Journal of Computational Physics 421 (2020) 109736.
  * [7] K. Carlberg, Y. Choi, S. Sargsyan, Conservative model reduction for finite-volume models, Journal of Computational Physics 371 (2018) 280–314.
  * [8] A. Schein, K. T. Carlberg, M. J. Zahr, Preserving general physical properties in model reduction of dynamical systems via constrained-optimization projection, International Journal for Numerical Methods in Engineering 122 (14) (2021) 3368–3399.
  * [9] K. Carlberg, M. Barone, H. Antil, Galerkin v. least-squares Petrov–Galerkin projection in nonlinear model reduction, Journal of Computational Physics 330 (2017) 693–734.
  * [10] M. Xiao, P. Breitkopf, R. F. Coelho, C. Knopf-Lenoir, P. Villon, W. Zhang, Constrained proper orthogonal decomposition based on QR-factorization for aerodynamical shape optimization, Applied Mathematics and Computation 223 (2013) 254–263.
  * [11] B. Sanderse, Energy-conserving discretization methods for the incompressible Navier-Stokes equations: application to the simulation of wind-turbine wakes, Ph.D. thesis, Eindhoven University of Technology (2013).


# Development of an Estimation Method for the Seismic Motion Reproducibility
of a Three-dimensional Ground Structure Model by combining Surface-observed
Seismic Motion and Three-dimensional Seismic Motion Analysis

Tsuyoshi Ichimura  Earthquake Research Institute and Department of Civil
Engineering, The University of Tokyo, Japan  Kohei Fujita  Earthquake Research
Institute and Department of Civil Engineering, The University of Tokyo, Japan
Ryota Kusakabe  Earthquake Research Institute and Department of Civil
Engineering, The University of Tokyo, Japan  Hiroyuki Fujiwara  Multi-hazard
Risk Assessment Research Division, National Research Institute for Earth
Science and Disaster Resilience, Japan  Muneo Hori  Research Institute for
Value-Added-Information Generation, Japan Agency for Marine-Earth Science and
Technology, Japan  Maddegedara Lalith  Earthquake Research Institute and
Department of Civil Engineering, The University of Tokyo, Japan

###### Abstract

The ground structure can substantially influence seismic ground motion
underscoring the need to develop a ground structure model with sufficient
reliability in terms of ground motion estimation for earthquake damage
mitigation. While many methods for generating ground structure models have
been proposed and used in practice, there remains room for enhancing their
reliability. In this study, amid many candidate 3D ground structure models
generated from geotechnical engineering knowledge, we propose a method for
selecting a credible 3D ground structure model capable of reproducing observed
earthquake ground motion, utilizing seismic ground motion data solely observed
at the ground surface and employing 3D seismic ground motion analysis. Through
a numerical experiment, we illustrate the efficacy of this approach. By
conducting $10^{2}$–$10^{3}$ cases of fast 3D seismic wave propagation
analyses using graphic processing units (GPUs), we demonstrate that a credible
3D ground structure model is selected according to the quantity of seismic
motion information. We show the effectiveness of the proposed method by
showing that the accuracy of seismic motions using ground structure models
that were selected from the pool of candidate models is higher than that using
ground structure models that were not selected from the pool of candidate
models.

##  1 Introduction

Structures on or within the ground can experience substantial impacts from the
ground structure during seismic events [1][2]. For instance, soft ground atop
the hard ground, coupled with the characteristics of seismic motion and the
ground structure, can lead to localized amplification of seismic motion,
resulting in damage. Enhancing the reliability of ground structure assessments
is required to enable accurate evaluation of structural behavior during
seismic events and mitigate earthquake-induced damage effectively.

Many methods for estimating ground structures have been proposed and employed
in practical applications. For instance, methods directly assessing in situ
geotechnical properties, such as borehole testing, yield reliable outcomes.
Nonetheless, due to the limited number of measurement points, interpolation
becomes necessary to estimate the spatial extent of the ground structure. In
the construction of 3D ground structure models, interpolation methods based on
geotechnical engineering knowledge, such as inverse distance weighting methods
[3], curvature minimization principles [4], and Kriging methods [5], are
employed. These interpolations rely on various assumptions, resulting in
multiple 3D ground structure models, necessitating additional evaluation to
determine the reliability of each model.

Such evaluation can utilize the observed ground motion of small amplitudes.
For instance, the reliability of a candidate ground structure model can be
assessed by comparing the phase velocity derived from long-term microtremor
observations [6] or by comparing the Green’s functions obtained via seismic
interferometry [7]. Alternatively, the quality of 3D ground structure models
can be evaluated by directly utilizing observed seismic ground motions of
small amplitudes, which typically exhibit stronger signal strengths. For
example, [8] demonstrates the feasibility of evaluating a 3D ground structure
model using observed seismic ground motions of small amplitudes. However, the
method in [8] assumes that seismic motions input to the 3D ground structure
model are observed at underground observation points, making it challenging to
apply to sites where underground observation points cannot be installed.

Building upon the above insights, we propose a method that selects a credible
ground structure model from many generated 3D ground structure models using
small amplitude seismic motions solely observed at the ground surface. We
illustrate the effectiveness of the method through numerical experiments by
conducting $10^{2}$–$10^{3}$ cases of fast 3D seismic wave propagation
analysis on GPUs, and show that we can select a credible ground structure
model based on the amount of seismic motion information observed at the ground
surface. Additionally, we show that the selected 3D ground structure model can
be utilized to evaluate ground motion with sufficient accuracy.

##  2 Method

Figure 1: Target system (left) and its numerical analysis model (right).

We introduce a method for extracting a credible 3D ground structure model from
a pool of candidate 3D ground structure models using small amplitude seismic
motions observed at the ground surface. The general setting, depicted on the
left of Fig. 1, involves inputting seismic waves to a 3D ground structure
comprising soft and hard soil, with resulting seismic ground motions observed
at the designated points on the surface. Notably, the 3D ground structure and
the input seismic waves are typically unknown.

Next, we elucidate the numerical analysis model for the target system, as
depicted on the right side of Fig. 1. This study uses the cartesian coordinate
system $x_{1},x_{2},x_{3}$ for simplicity. The ground enclosed by the dashed
line in the left panel of Fig. 1 is designated as $V(\mathbf{m})$, defined by
a set of model parameters $\mathbf{m}$. An external force, represented by the
input velocity wave $f_{i}(t)$, is applied to the bottom of $V(\mathbf{m})$,
and its dynamic response is analyzed to yield the time-history response
$obs^{k}_{i}(t)$ at each observation point $k$ on the ground surface. To
compute the small amplitude seismic response of the ground to an earthquake,
we model $V(\mathbf{m})$ as linearly elastic and solve the governing equations
for a linear dynamic elastic body

$d_{i}(c_{ijkl},d_{k}u_{l})=\rho\ddot{u}_{j}.$ (1)

Here, $c_{ijkl}$, $u_{j}$, and $\rho$ represent the elasticity tensor,
displacement in the $j$-th direction, and density, respectively. The notation
$(\ddot{~{}})$ denotes second-order derivatives in time, and $d_{i}$ indicates
the differentiation in the $i$-th direction. Note that stress-free, non-
reflective, and semi-infinite absorbing boundary conditions are enforced on
the model’s top, bottom, and side surfaces. We aim to develop a method that
can evaluate $V(\mathbf{m})$ and $f_{i}(t)$ that is consistent with
$obs^{k}_{i}(t)$.

We describe the parametrization employed for this evaluation. First, we
suppose that the input velocity wave can be described as follows:

$f_{i}(t)=\sum_{j}c_{ij}p(t-(j-1)\Delta t).$

We utilize unit impulse waves $p(t)$, for instance, trigonometric functions,
to represent the input velocity wave. Here, $c_{ij}$ represents unknown scalar
values assessed through error minimization, and $\Delta t$ denotes the time-
stepping stride for discretization. While arbitrary $p(t)$ and $\Delta t$ are
permissible, we configure them to facilitate the reconstruction of the input
wave within the target frequency range. Moreover, we denote the Green’s
function of the $i$-th directional velocity response at observation point $k$
when $p(t)$ is input in the $j$-th direction at the bottom of $V(\mathbf{m})$,
as $G^{k}_{ij}(\mathbf{m},t)$. The $i$-th directional velocity at observation
point $k$ for $f_{j}(t)$ can thus be expressed as

$U^{k}_{i}(\mathbf{m},t)=\sum_{j=1}^{3}\sum_{l}G_{ij}^{k}(\mathbf{m},t-(l-1)\Delta
t)c_{jl}.$

Utilizing the observed ground velocity in the $i$-th direction at the $k$-th
observation point ($obs^{k}_{i}(t)$), the error can be assessed as follows:

$ERR(\mathbf{m})=\frac{1}{3n_{obs}}\sum_{k=1}^{n_{obs}}\sum_{i=1}^{3}\frac{\sqrt{\int(U^{k}_{i}(\mathbf{m},t)-obs^{k}_{i}(t))^{2}\mathrm{d}t}}{\sqrt{\int
obs^{k}_{i}(t)^{2}\mathrm{d}t}}.$ (2)

With the parameter settings outlined above, we compute $c_{ij}$ that minimizes
the error from the observed seismic motion $obs^{k}_{i}(t)$ for each candidate
$l$-th 3D ground structure model characterized by model parameters
$\mathbf{m}_{l}$. First, we compute $G^{k}_{ij}(\mathbf{m}_{l},t)$ for the
model parameters $\mathbf{m}_{l}$. Subsequently, by considering the stationary
condition of $ERR$, which is a second-order function of $c_{ij}$, the
coefficients $c_{ij}$ can be determined by solving:

$\mathbf{A}\mathbf{c}=\mathbf{b}.$

Here, $\mathbf{A}$, $\mathbf{b}$ represent a constant matrix and a constant
vector, respectively, while $\mathbf{c}$ denotes an unknown vector with
components $c_{ij}$.

By applying singular value decomposition to $\mathbf{A}$ and discarding
negligible singular values, we construct a pseudoinverse matrix and compute
$\mathbf{c}$ and $ERR$. Note that the dimension of $\mathbf{A}$ is small
because this computation can be performed independently for each trial, and
thus the cost of this computation can be kept small enough. That is, if the
seismic ground motions contain sufficient information (i.e., the number of
observation points and the number of events is adequate) and are suitably
constrained by the model-specific time-history Green’s functions, it is
possible to assess the credibility of various 3D ground structure models by
comparing the error $ERR$ obtained through attempts to reproduce observed
seismic ground motions using $\mathbf{m}_{l}$.

##  3 Numerical Experiment

We conducted a numerical experiment to demonstrate the efficacy of the
proposed method in selecting a credible 3D ground structure model from among
geotechnically estimated ground models, using seismic ground motions observed
at the ground surface. Here, a credible ground model can reproduce the seismic
ground motions observed at the surface. Specifically, the method described in
Section 2 estimates the seismic ground motion that best fits the observed
seismic ground motion using a given 3D ground structure model, albeit with
introduced errors. This analysis in estimating seismic ground motion that best
fits the observed ground motion for a given 3D ground structure model is
conducted each time an earthquake is observed. Discarding ground models that
inadequately reproduce seismic ground motions can identify a credible ground
model capable of reproducing observed seismic ground motions at a specific
site.

First, we set a reference model to serve as the correct solution for the
numerical experiment. The reference model comprises a two-layer structure
comprising a sedimentary layer and bedrock, which imitates a real ground
structure. It spans 600, 600, and 100 m in the $x_{1}$, $x_{2}$, and $x_{3}$
directions, respectively (see Fig. 2). The ground surface is flat at an
elevation of 0 m, with the thickness of the first layer specified as depicted
in Fig. 7a). The physical properties of each layer are given in Fig. 2. This
ground structure emulates valley floor lowlands formed by sediment
accumulation from river meandering. During earthquakes, local ground
amplification occurs, emphasizing the importance of accurately assessing the
ground structure for evaluating seismic ground motion and making effective
earthquake mitigation strategies. Assuming that the reference model is
situated at KiK-net [9] station IBRH19, we simulate a scenario where
earthquakes listed in Table 1 occur sequentially, and seismic motions are
observed. To simulate this scenario, actual ground motions observed at
underground stations of IBRH19 were inputted as forces from the bottom of the
reference ground model, and pseudo-observed seismic ground motions were
recorded at the ground surface. The analysis was performed up to 2.5 Hz, which
is the frequency range considered to have a significant impact on structural
damage. Semi-infinite absorbing boundary conditions are applied to the sides
and bottom of the numerical model.

Figure 2: Reference ground model. The model comprises two layers, with a flat surface and a sedimentary layer with varying thickness. The thickness of the sedimentary layer is illustrated in Fig. 7a). Table 1: Properties of the earthquakes used in the numerical experiment Event # | Data time (JST) | Epicenter | Depth (km) | Magnitude
---|---|---|---|---
1 | 2023/02/25 22:27 | 42.755N 145.075E | 63 | 6.0
2 | 2023/03/27 00:04 | 38.307N 141.615E | 60 | 5.3
3 | 2023/05/11 04:16 | 35.170N 140.185E | 40 | 5.2
4 | 2023/07/29 19:34 | 36.347N 139.958E | 77 | 4.6
5 | 2023/12/22 01:11 | 35.238N 141.137E | 10 | 4.8
6 | 2024/01/28 08:59 | 35.6N 140.0E | 80 | 4.8

Next, a set of candidate ground structure models is generated using
geotechnical engineering methods. While many methods generating ground
structure models exist, we assume that the physical properties of the first
and second layers, along with the boundary location between them, are
determined from a simple ground survey conducted at the 120 survey points
illustrated in Fig. 7a) (i.e., $\mathbf{m}$ represents model parameters
describing the boundary shape between the first and second layers). The
geometry of the boundary between the first and second layers within the target
area of 600 $\times$ 600 m is estimated based on information regarding the
boundary layer location at these 120 survey points, creating a set of
candidate ground structure models. Although many methods for estimating the
layer boundary geometry under these conditions exist, we employ an inverse
distance weighting method. This method assigns weights to observed values
obtained near the evaluation point based on the inverse of the distance,
yielding an interpolated result. Utilizing the inverse distance weighting
method, the layer thickness $z$ at any point in the domain is evaluated as
follows:

$z(x_{1},x_{2})=\sum_{i=1}^{M}\frac{\bar{z}_{i}/d_{i}^{q}}{\sum_{i=1}^{M}1/d_{i}^{q}},~{}{\rm
where}~{}d_{i}=\sqrt{(x_{1}-\bar{x}_{1}^{i})^{2}+(x_{2}-\bar{x}_{2}^{i})^{2}},$

where $\bar{z}_{i}$ denotes the layer thickness at the $i$-th measurement
point located at $(\bar{x}_{1}^{i},\bar{x}_{2}^{i})$, $M$ represents the
number of measurement points nearest to the target position $(x_{1},x_{2})$
utilized for interpolation, and $q$ is a positive constant parameter. The
values of $M$ and $q$ are arbitrary. In this context, we designate
$M=i~{}(i=1,2,...,50)$ and $q=0.1i~{}(i=1,2,....,40)$ to generate a total of
$50\times 40=2000$ 3D ground structure models. Subsequently, from the pool of
candidate 3D ground structure models, we select 236 models that exhibit
substantial variations in the layer boundary geometry, based on a 3D ground
structure model generated using commonly used parameters (i.e., $M=20$ and
$q=2.0$). Note that Laplace smoothing was applied five times to mitigate steep
slopes in certain areas.

Utilizing the method described in Section 2 and the pseudo-observed seismic
ground motions described above, we attempt to identify credible 3D ground
structure models by reproducing the observed seismic ground motions for each
of the 236 obtained models. First, we consider a scenario where only one
observation point is available on the ground surface at
$(x_{1},x_{2})=(300,300)$ m. The analysis conducted for each candidate model
is the same as the analysis using the reference model (i.e., the target
frequency is set up to 2.5 Hz), with semi-infinite absorbing boundary
conditions applied to the sides and bottom of the model (the same applies to
the other analysis cases in this Section). Figure 3a) illustrates the errors
obtained for each ground structure model across each earthquake event.
Although the error varies depending on the model, it remains very small
regardless of the ground structure model in the case of a single observation
point (from the definition of $ERR$ in Eq. (2), $ERR=0.05$ means that the
waveform is matched within an error of about 5% of its amplitude). This
indicates that the waveform constraints are insufficient when the number of
observation points is limited. In other words, despite the complexity of the
time history Green’s function in a 3D medium, observed waveforms can be
accurately reconstructed regardless of the model used by imposing errors on
the estimated input waveform. Indeed, the estimated input waveform obtained
using model000083 (i.e., model number 83 of the 236 candidate models), which
exhibits a small error, substantially differs from the true waveform in
amplitude and phase characteristics (Fig. 4). Conversely, increasing the
number of observation points becomes imperative to leverage model-specific
time-history constraints in Green’s functions for accurately reconstructing
both observed and input waveforms. This constraint can be leveraged to
diminish the ability to reproduce the observed waveforms depending on the
ground structure model, thereby facilitating the selection of the appropriate
ground structure model.

a) $ERR$ for each event (left: overview, right: closeup)

b) Maximum $ERR$ for events (left: overview, right: closeup)

Figure 3: Estimated $ERR$ using one observation point

a) True incident wave

b) Estimated incident wave using model000083 ($ERR=0.0243$)

Figure 4: Estimated incident wave using one observation point for event #1.
Estimation accuracy is low even if $ERR$ is low.

Next we employ nine observation points on the ground surface
$(x_{1},x_{2})=(100+200i,100+200j)$ m, where $i,j=0,1,2$. Figure 5a)
illustrates the error obtained for each model across each earthquake event.
Compared with the single observation point case, the information derived from
the pseudo-observed earthquake ground motion adequately constrains the ground
model. Consequently, models exhibiting small error levels are consistently
identified irrespective of the earthquake event. The comparison with the
historical maximum error depicted in Fig. 5b) reveals that the credible ground
model remains consistently selected even as the number of experienced
earthquakes and the quantity of information increases. Furthermore, we augment
the information by increasing the number of surface observation points to 25,
located at $(x_{1},x_{2})=(100+100i,100+100j)$ m, where $i,j=0,1,....,4$. The
error estimation is shown in Fig. 6a). As in the 9-point scenario, models with
small errors are systematically chosen. The comparison with the historical
maximum error exhibited in Fig. 6b) underscores the continued stable selection
of the credible ground model even with the heightened number of experienced
earthquakes and augmented information. Comparing the 9-point case with the
25-point case, it is evident that the increase in information facilitates the
stable and systematic extraction of ground models characterized by smaller
errors.

a) $ERR$ for each event (left: overview, right: closeup)

b) Maximum $ERR$ for events (left: overview, right: closeup)

Figure 5: Estimated $ERR$ using nine observation points

a) $ERR$ for each event (left: overview, right: closeup)

b) Maximum $ERR$ for events (left: overview, right: closeup)

Figure 6: Estimated $ERR$ using 25 observation points

Finally, we checked the extracted results. Figures 7a,b, and c) display the
layer thickness distribution of the reference ground model, model000018 (with
a minimum error of 0.388 for the case with 25 observation points), and
model000228 (with error of 0.781 for the case with 25 observation points).
Model000018, with the smallest error, closely resembles the reference ground
model, indicating the extraction of a plausible 3D ground structure model
through this method. While no model among the candidate ground models exactly
matches the reference model, model000018 was selected due to its minimal
difference compared to the target wavelength used for analysis. Conversely,
model000228, with a large error, substantially deviates from the reference
model. Next, Fig. 8 illustrates the $x_{1}$ component of the input seismic
motion estimated by model000018 and model000228 using pseudo-observed seismic
ground motion at 25 observation points during earthquake event #1. The input
waveform estimated using model000018, evaluated to have a small error, closely
resembles the true input wave, while the waveform estimated using model000228
differs substantially from the true input wave. This underscores that
sufficient observed seismic ground motions enable the Green’s function
constraint of the time history in the 3D medium obtained from a 3D wavefield
analysis, ensuring estimation performance for both observed and input seismic
motions when using a plausible model. Finally, Fig. 9 depicts the ground
motion distribution using the reference model, model000018, and model000228.
The responses of the reference model and model000018 closely align regardless
of the location from the observation points, whereas the response of
model000228 differs markedly. As shown in Fig. 10, model000018, selected as
the most credible 3D ground structure model in this study, can estimate not
only the time series response at the observation point
($(x_{1},x_{2})=(300,300)$ m) but also the time series response at a point
which is located away from the observation points ($(x_{1},x_{2})=(350,350)$
m). This indicates that a consistent ground model that matches the observed
seismic ground motion has been extracted. These results underscore the
significance of selecting a plausible 3D ground structure model, as
demonstrated in this study, for mitigating earthquake damage, since even 3D
ground structure models estimated by geotechnical engineering methods exhibit
substantial differences in their performance in reproducing seismic ground
motions.

a) Reference model

b) Model000018

c) Model000228

Figure 7: Thickness of sedimentary layers for each model. The position of 120
boring log points is also indicated in a).

As described above, it is evident that a 3D ground structure model
demonstrating high performance in reproducing observed seismic motions can be
extracted by using seismic motion observations at the ground surface.
Achieving this necessitates multiple 3D seismic response analyses. In this
study, finite element models are generated for the 236 ground models with
varying layer thicknesses, and for each model, the Green’s function response
of Eq. (1) to the unit impulse wave $p(t)$ in each of the $x_{1}$, $x_{2}$,
and $x_{3}$ directions are computed, resulting in a total of 236$\times$3 =
708 instances of 3D wave field analysis. An automatic method for generating
high-quality finite element models is required to facilitate such a large
number of analyses. Moreover, a fast implicit time integration-based 3D
seismic analysis method capable of performing stable calculations, even on
finite element models with small local elements that are essential to
faithfully model the complex ground geometry, is required. Consequently, in
this study, finite element models comprising second-order tetrahedral elements
are automatically generated using a mesh generator [10] employing octree
background cells (the finite element models utilized in this study average
2,320,401 degrees of freedom). The Newmark-$\beta$ method, a type of implicit
time integration, is employed to compute the seismic ground motion over 16,000
time steps with a time increment of $dt=0.001$ s to ensure stability. While
implicit time integration is suitable for stable computation, it entails
solving solutions of large sparse matrix equations, resulting in substantial
computational costs. To address this issue, we use a fast solver from [11],
based on the conjugate gradient method with variable preconditioning [12].
Consequently, even a 3D seismic analysis with 2.32 million degrees of freedom
(516748 second-order tetrahedral elements with minimum element size of 5 m)
and 16,000 time steps can be executed in approximately 270 s on a computing
environment with an NVIDIA A100 PCIe 40 GB GPU [13]. The entire evaluation can
be completed stably within a short duration of approximately 1.5 h using 44
A100 GPUs. Due to time constraints, conducting 3D wave propagation analysis
for numerous ground models of this scale would have been impractical in the
past. However, leveraging a fast 3D wave propagation method, as demonstrated
here, enables such analyses to be conducted within feasible timeframes.

a) True incident wave

b) Estimated incident wave using model000018 ($ERR=0.388$)

c) Estimated incident wave using model000228 ($ERR=0.781$)

Figure 8: Estimated incident wave using 25 observation points for event #1.
The incident wave can be estimated accurately by selecting a model with a low
$ERR$.
Figure 9: Displacement response at the surface for the true incident wave
(event #1) at $t=50$ s

a) Wave at an observation point at $(x_{1},x_{2})=(300,300)$ m.

b) Wave at a point on the surface at $(x_{1},x_{2})=(350,350)$ m, which is far
from the observation points.

Figure 10: Displacement time-history response for the true incident wave
(event #1) using the reference model (black) and model000018 (blue).

##  4 Concluding Remarks

This study introduced a method to extract a 3D ground structure model that can
reproduce observed seismic ground motions from a pool of candidate ground
structure models using surface-observed earthquake ground motions and 3D
seismic ground motion analysis. Through numerical experiments, we illustrated
its effectiveness. Even when ground models are generated using methods based
on geotechnical engineering aspects, their performance in reproducing observed
ground motion varies substantially. A comparison of simulation results between
models exhibiting high and low performance underscored substantial differences
in seismic damage estimation during earthquakes, highlighting the efficacy of
our method. While our numerical experiments focused solely on treating the
layer boundary geometry as a model parameter, there are no constraints on
parameterization in this context. Both the material properties and geometry
can serve as model parameters for extracting a plausible 3D ground structure
model. The above findings indicate that our proposed method is expected to
contribute to enhancing the reliability of 3D ground structure models
generated by supplementing borehole data and other databases, which is
expected to become more common in the future. There is also possibility in
enhancing the method by utilizing further HPC and AI for efficiently
generating better candidate models around the ground models that are chosen to
be credible and searching for a ground model with better performance in
reproducing observed ground motions.

### Acknowledgments

This work used waveforms provided by NIED KiK-net, National Research Institute
for Earth Science and Disaster Resilience, Japan.

## References

  * [1] Liang, J., Sun S.: Site Effects on Seismic Behavior of Pipelines: A Review, ASME J. Pressure Vessel Technol., 122, pp. 469-475 (2000)
  * [2] Advisory Committee for Natural Resources and Energy, Urban Area Thermal Energy Committee, Gas Safety Subcommittee, Working Group for Earthquake Disaster Prevention: Report on Disaster Mitigation for Gas Supply in View of Great East Japan Earthquake (in Japanese) (2012)
  * [3] Shepard, D.: A two-dimensional interpolation function for irregularly-spaced data. In Proceedings of the 1968 23rd ACM national conference (ACM ’68). Association for Computing Machinery, New York, NY, USA, 517–524. https://doi.org/10.1145/800186.810616 (1968)
  * [4] Briggs, I. C.: Machine contouring using minimum curvature. Geophysics; 39 (1): 39–48. doi: https://doi.org/10.1190/1.1440410 (1974)
  * [5] Oliver, M. A., Webster, R.: Kriging: a method of interpolation for geographical information systems, International journal of geographical information systems, 4:3, 313-332, DOI: 10.1080/02693799008941549 (1990)
  * [6] Aki, K., Space and time spectra of stationary stochastic waves with special reference to microtremors. Bull. Erthq. Res. Inst, 35, 415–456 (1957)
  * [7] Wapenaar, K., Fokkema, J.: Green’s function representation for seismic interferometry, Geophysics, Vol.71, No.4, pp. SI33–SI46 (2006)
  * [8] Yamaguchi, T., Ichimura, T., Fujita, K., Hori, M., Wijerathne, L., Ueda, N.: Data-Driven Approach to Inversion Analysis of Three-Dimensional Inner Soil Structure via Wave Propagation Analysis. In: Krzhizhanovskaya, V.V., et al. Computational Science - ICCS 2020. ICCS 2020. Lecture Notes in Computer Science, 12139. Springer, Cham (2020)
  * [9] National Research Institute for Earth Science and Disaster Resilience: NIED K-NET, KiK-net, National Research Institute for Earth Science and Disaster Resilience, doi:10.17598/NIED.0004 (2019)
  * [10] Ichimura, T., Hori, M., and Bielak, J.: A hybrid multiresolution meshing technique for finite element three-dimensional earthquake ground motion modelling in basins including topography, Geophysical Journal International, 177(3), 1221–1232 (2009)
  * [11] Kusakabe, R., Fujita, K., Ichimura, T., Yamaguchi, T., Hori, M., Wijerathne, L.: Development of regional simulation of seismic ground-motion and induced liquefaction enhanced by GPU computing, Earthquake Engineering & Structural Dynamics, 50, 197–213 (2021)
  * [12] Ichimura, T., Fujita, K., Tanaka, S., Hori, M., Maddegedara, L., Shizawa, Y., Kobayashi, H.: Physics-Based Urban Earthquake Simulation Enhanced by 10.7 BlnDOF $\times$ 30 K Time-Step Unstructured FE Non-Linear Seismic Wave Simulation, SC ’14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, New Orleans, LA, USA, 2014, pp. 15–26 (2014)
  * [13] NVIDIA Ampere architecture whitepaper, https://www.nvidia.com/content/dam/en-zz/Solutions/Data-Center/nvidia-ampere-architecture-whitepaper.pdf, last accessed 2024/2/27