Bathymetry & Boundary Conditions
================================
In this part of our project, we will first introduce bathymetry to our one-dimensional solver.
Then, we will simulate a one-dimensional flow over a weir and study hydraulic jumps.
Last, we run a one-dimensional tsunami simulation.
.. _ch:bathymetry:
Non-zero Source Term
--------------------
Bathymetry is the topography of the ocean and the most important part of the shallow water equation's source term when simulating tsunamis.
We use sea level as reference: Bathymetry in a dry cell (onshore) :math:`\mathcal{C}_i` is zero or positive: :math:`b_i \ge 0`, a wet cell :math:`\mathcal{C}_j` (offshore) has a negative value: :math:`b_j < 0`.
The introduction of bathymetry into the f-wave solver is comparably simple because it can be incorporated directly into the decomposition of the jump in the flux function:
.. math:: \Delta f - \Delta x \Psi_{i-1/2} = \sum_{p=1}^2 Z_p.
:label: eq:fwave_source
The term :math:`\Delta x \Psi_{i-1/2}` summarizes\ [1]_ the effect of the bathymetry:
.. math::
:label: eq:psi
\Delta x \Psi_{i-1/2} := \begin{bmatrix}
0 \\
-g (b_r - b_l) \frac{h_l+h_r}{2}
\end{bmatrix}.
.. _ch:reflecting_boundary_conditions:
.. admonition:: Tasks
#. Extend the f-wave solver with support for bathymetry according to :eq:`eq:fwave_source` and :eq:`eq:psi`.
#. Implement an example which illustrates the effect of bathymetry in your extended solver.
**Note:** You have to add bathymetry to the respective patches implementing the abstract class ``patches::WavePropagation`` as well.
Also, ghost cells require a valid value for the bathymetry, i.e., set :math:`b_0 := b_1` and :math:`b_{n+1} = b_n`.
.. [1]
Details can be found in: *A wave propagation method for conservation
laws and balance laws with spatially varying flux functions*. D. S.
Bale et. al., 2003.
Reflecting Boundary Conditions
------------------------------
The introduction of bathymetry naturally leads to the question of handling wetting and drying in near-shore areas, i.e., :math:`b(x) \approx 0`.
The f-wave solver is not capable of handling wetting and drying. We will therefore impose two simplifications to keep our solver stable:
* Cells which are initially wet stay wet for the entire simulation:
We assure each cell has enough water initially for the entire duration of a simulation.
A fix for tsunami simulations will be introduced in :numref:`ch:tsunami_1d`.
* Cells which are initially dry stay dry for the entire simulation:
At a wet-dry interface reflecting boundary conditions are used.
If a wave moves from a wet region towards the shore, we assume an infinitely high solid wall at the shore-line.
This is equivalent to reflecting boundary conditions.
Since no water is able to pass the shoreline, we have to ensure that the particle velocity at the boundary is zero.
.. figure:: data_2/wet_dry_boundary.svg
:name: fig:wet_dry_boundary
:align: center
:width: 35.0%
Wet-dry boundary for a wet cell :math:`\mathcal{C}_{i-1}` and dry
cell :math:`\mathcal{C}_i`.
Our finite volume scheme can have neighboring cells :math:`\mathcal{C}_{i-1}` and :math:`\mathcal{C}_i`, where :math:`\mathcal{C}_{i-1}` is wet and :math:`\mathcal{C}_{i}` is dry.
:numref:`fig:wet_dry_boundary` illustrates this setting.
As seen in the shock-shock example of :numref:`ch:shock_rare`, we obtain a middle state with zero particle velocity if we set an equal water height and an opposite-sign momentum in the dry cell.
Additionally, setting the bathymetry in the dry cell to the value of the wet cell completes the setup of reflecting boundary conditions:
.. math::
:label: eq:reflecting_bc
h_{i} &:= h_{i-1} \\
(hu)_{i} &:= -(hu)_{i-1} \\
b_{i} &:= b_{i-1}
.. admonition:: Tasks
#. Implement reflecting boundary conditions as defined in :eq:`eq:reflecting_bc`.
#. Show that you obtain the :numref:`ch:shock_rare`’s one-sided solution of the shock-shock setup if you set :math:`q_l` everywhere initially and use reflecting boundary conditions at the right boundary, and outflow boundary conditions at the left boundary.
**Note:** Reflecting boundary conditions for the right boundary can be implemented through a ghost cell :math:`\mathcal{C}_{n+1}` that is dry.
.. _ch:hydraulic_jumps:
Hydraulic Jumps
---------------
This section illustrates the behavior of flows in a sub- and supercritical example.
A lab-experiment is recorded in the video `Hydraulics over a weir <https://serc.carleton.edu/details/files/19076.html>`_.
In shallow water theory we can characterize this behavior by the Froude number:
.. math:: F := \frac{u}{\sqrt{gh}}.
We call regions with :math:`F < 1` subcritical, :math:`F \approx 1` critical and :math:`F > 1` supercritical.
This definition matches the wave speeds of :numref:`ch:riemann_solver` in terms of the eigenvalues :math:`\lambda_{1/2} = u \mp \sqrt{gh}`.
In the supercritical case the wave velocity relative to the fluid, i.e., :math:`\mp\sqrt{gh}`, is smaller than the velocity of the flow.
Thus, information travels in one direction only and hydraulic jumps can occur.
This behavior of flow is closely related to the subsonic, sonic and supersonic cases in acoustics.
The two following examples\ [2]_ set up a sub- and supercritical flow over a hump.
We use the interval :math:`(0, 25)` as computational domain and outflow boundary conditions everywhere.
For the **subcritical flow** we use the following initial values:
.. math::
:label: eq:subcritical_example1
\begin{aligned}
b(x) &=
\begin{cases}
-1.8 - 0.05 (x-10)^2 \quad &\text{if } x \in (8,12) \\
-2 \quad &\text{else}
\end{cases}\\
h(x, 0) &= -b(x) \quad \text{if } x \in [0,25] \\
hu(x, 0) &= 4.42 \quad \text{if } x \in [0,25].
\end{aligned}
For the **supercritical flow** we use the following initial values:
.. math::
:label: eq:supercritical_example1
\begin{aligned}
b(x) &=
\begin{cases}
-0.13 - 0.05 (x-10)^2 \quad &\text{if } x \in (8,12) \\
-0.33 \quad &\text{else}
\end{cases}\\
h(x, 0) &= -b(x) \quad \text{if } x \in [0,25] \\
hu(x, 0) &= 0.18 \quad \text{if } x \in [0,25].
\end{aligned}
.. _tasks-2-1:
.. admonition:: Tasks
#. Compute the location and value of the maximum Froude number for the subcritical setting given in :eq:`eq:subcritical_example1` and the supercritical setting given in :eq:`eq:supercritical_example1` at the initial time :math:`t=0`.
#. Implement both cases through the base class ``setup::Setup.h``.
:math:`t \in [0, 200]` is a reasonable time window for your simulation.
#. Determine the position of the hydraulic jump (stationary discontinuity) in your supercritical solution and show that our f-wave solver fails to converge to the analytically expected constant momentum over the entire domain.
.. [2]
Details of the original benchmarks can be found in: *Augmented
Riemann solvers for the shallow water equations over variable
topography with steady states and inundation*. D. L. George, 2008.
.. _ch:tsunami_1d:
1D Tsunami Simulation
---------------------
Now, with support for bathymetry in our Riemann solver, we can run our first one-dimensional simulations using field data.
For this, we first extract input data for our solver from a Digital Elevation Model (DEM).
.. figure:: data_2/1d_domain.png
:name: fig:1d_domain
:align: center
:width: 100%
Visualization of the GEBCO_2021 Grid and the one-dimensional domain as a line between the two points :math:`p_1=(141.024949, 37.316569)` and :math:`p_2=(146.0, 37.316569)`. The epicenter of the March 11, 2011, M 9.1 Tohoku earthquake is north of the considered domain.
As illustrated in :numref:`fig:1d_domain`, we use the `GEBCO_2021 Grid <https://www.gebco.net/data_and_products/gridded_bathymetry_data/>`_ to extract the bathymetry profile between :math:`p_1=(141.024949, 37.316569)` and :math:`p_2=(146.0, 37.316569)`.
The `Generic Mapping Tools <https://www.generic-mapping-tools.org/>`_' function `grdtrack <https://docs.generic-mapping-tools.org/latest/grdtrack.html>`_ is able to extract a 1D profile from a DEM.
:numref:`fig:1d_profile` shows the bathymetry data required by our solver for the one-dimensional domain.
.. figure:: data_2/1d_profile.png
:name: fig:1d_profile
:align: center
:width: 100%
Visualization of the required bathymetry data, extracted from the GEBCO_2021 Grid.
Next, we extend our software by adding a simple CSV reader, which can parse the data for us.
This reader is then used in a new setup ``setups::TsunamiEvent1d``.
The setup sets the following initial values for the water height :math:`h`, momentum :math:`hu`, and bathymetry :math:`b`:
.. math::
:label: eq:tsunami_event_1d
\begin{split}
h &= \begin{cases}
\max( -b_\text{in}, \delta), &\text{if } b_\text{in} < 0 \\
0, &\text{else}
\end{cases}\\
hu &= 0\\
b &= \begin{cases}
\min(b_\text{in}, -\delta) + d, & \text{ if } b_\text{in} < 0\\
\max(b_\text{in}, \delta) + d, & \text{ else}.
\end{cases}
\end{split}
:math:`b_\text{in}(x) \in \mathbb{R}` is the bathymetry extracted from the DEM.
:math:`d(x) \in \mathbb{R}` is the vertical displacement, typically caused by a subduction-zone earthquake.
In this task, we'll use artificial displacement data, i.e., we assume:
.. math::
d(x) = \begin{cases}
10\cdot\sin(\frac{x-175000}{37500} \pi + \pi), & \text{ if } 175000 < x < 250000 \\
0, &\text{else}.
\end{cases}
The constant :math:`\delta \in \mathbb{R}^+` avoids running into numerical issues due to missing support for wetting and drying in our solver.
We will typically use :math:`\delta := 20\,\text{m}` and apply reflecting boundary conditions whenever a wet cell neighbors a dry cell.
.. admonition:: Tasks
#. Extract bathymetry data for the 1D domain, shown in :numref:`fig:1d_domain`, from the GEBCO_2021 Grid. Use a 250m sampling between the two points :math:`p_1` and :math:`p_2`.
#. Extend the class ``tsunami_lab::io::Csv`` such that it can read your extracted bathymetry data.
#. Add a new setup ``setups::TsunamiEvent1d`` which uses your CSV reader and the artificial displacement to initialize your quantities as defined in :eq:`eq:tsunami_event_1d`.
#. Run and visualize the setup. What runup do you observe?
#. (optional) Modify the displacement! What are the impacts?