Settings and parameters

Kestrel input files are divided into different blocks, which may be specified in any order. We document the options for each below.

The settings are given in the form keyword = value. Here we provide example values for each keyword.

Some settings are required and the simulation will not commence if values are not given.

Other settings are optional and default values will be used if they are not given in the input file. In this case the default value is used as the example.

There are also some settings that are required in only some circumstances (conditionally required), and some that are optionally used in some circumstances (conditionally optional).

Domain

The Domain block specifies the location, size, resolution and boundary conditions that are applied at domain edges. The domain block is identified using the block keyword Domain:.

The domain is an interval in one-dimensional simulations and a rectangular region for two-dimensional simulations. In one dimension, the interval is taken to be along the \(x\) axis. In two dimensions, the coordinates are \(x\) and \(y\). For georeferenced simulations on topography, the domain is rectangular in the WGS84 UTM zone corresponding to the centre of the domain, with \(x\) corresponding to the Easting coordinate and \(y\) corresponding to the Northing coordinate.

The domain is divided into equal sized tiles, with the tiles divided into equal sized cells; this allows for efficient computation for flows following irregular paths along topography.

The following required settings are specified in the Domain block:

Lat = 13.248745

The latitude of the centre of the domain in decimal degrees in WGS84 coordinates, for georeferenced simulations. For simulations using an analytical surface, this should be set to the required central point for \(y\)-coordinate.

Lon = 123.686939

The longitude of the centre of the domain in decimal degrees in WGS84 coordinates, for georeferenced simulations. For simulations using an analytical surface, this should be set to the required central point for \(x\)-coordinate.

nXtiles = 40

Number of tiles in the \(x\)-direction.

nYtiles = 10

Number of tiles in the \(y\)-direction.

nXpertile = 25

Number of finite volume cells in the \(x\)-direction in each tile.

nYpertile = 25

Number of finite volume cells in the \(y\)-direction in each tile.

Xtilesize = 50

Length of a tile in the \(x\)-direction in metres. (Currently this also sets the length of a tile in the \(y\)-direction.)

Note

The numerical resolution of the simulation is determined by the cell size, given by Xtilesize / nXpertile.

Note

To perform one-dimensional simulations, set nYtiles = 1 and nYpertile = 1.

The following optional settings can be specified (here given with their default value if left unset):

Boundary conditions = halt

Sets the conditions that apply at the edge of the domain. Options are:

halt (default) – simulation terminates if flow reaches domain boundary;

sponge – flow quantities are forced to decay near to the domain boundary;

periodic – flow quantities at a boundary are transferred to opposite boundary;

Warning

Periodic boundary conditions do not enforce periodicity in the initial bed elevation. This means that applying them for non-periodic topographies (such as constant slopes) will lead to step changes in \(b\) across the boundaries.

dirichlet – flow quantities at a boundary are specified; values for \(H\), \(\bar{u}\), \(\bar{v}\) and \(\bar{\psi}\) are needed, and are specified using the domain block variables Boundary Hn, Boundary U, Boundary V, and Boundary psi, respectively.

The following settings are conditionally required when the optional Boundary conditions = dirichlet:

Boundary Hn

Value of flow depth, \(H\), to impose at a domain edge if Boundary conditions = dirichlet. Unused by default.

Boundary U

Value of flow velocity component in the \(x\)-direction, \(\bar{u}\), to impose at a domain edge if Boundary conditions = dirichlet. Unused by default.

Boundary V

Value of flow velocity component in the \(y\)-direction, \(\bar{v}\), to impose at a domain edge if Boundary conditions = dirichlet. Unused by default.

Boundary psi

Value of flow solids concentration, \(\bar{\psi}\), to impose at a domain edge if Boundary conditions = dirichlet. Unused by default.

Note

Dirichlet conditions only affect parts of the simulation containing flow. Therefore, the typical way to make use of them is to provide some flowing material via a Cap block (see below) that touches the edge of the domain.

Initial and source conditions

Flowing material may be input into simulations via Cap, Cube and Source blocks. Any number of these may be added to the input file.

Note

The Initial condition declaration, which is part of the Solver block (described below) may also be used to specify arbitrary initial conditions by providing valid solution data file(s).

Cap

A Cap block places a volume of flowing material with a circular base into the simulation at the initial time. The following declarations are required:

The location of the volume can be specified by giving either

• the latitude (capLat) and longitude (capLon) of its centre

or

• the offset of the volume’s centre from the centre of the domain (capX, capY), in metres.

These required specifiers give:

capLat = 13.248745

The latitude of the centre of the volume in decimal degrees in WGS84 coordinates.

capLat = 123.686939

The longitude of the centre of the volume in decimal degrees in WGS84 coordinates.

capX = 100

The offset of the centre of the volume along the \(x\) axis in metres from the centre of the domain.

capY = -50

The offset of the centre of the volume along the \(y\) axis in metres from the centre of the domain.

To set the volume of the release, exactly two out of three of the following settings are required:

capRadius

Sets the radius of the circular base in metres.

Note

The radius should be large enough to ensure that the volume can be represented on the numerical grid.

capHeight

Sets the height of the volume in metres. N.b. the precise meaning of this depends on the value of capShape (see below).

capVolume

Sets the volume in metres.

From the two specified values, Kestrel determines the third.

The following optional parameters set additional attributes for the volume release:

capU = 0

Sets the initial velocity \(\bar{u}\) of the release in the \(x\)-direction (m/s).

capV = 0

Sets the initial velocity \(\bar{v}\) of the release in the \(y\)-direction (m/s).

capConc = 0

Sets the volumetric fraction \(\bar{\psi}\) of solids present in the release.

capShape = para

This selects the geometry of the volume. There are three options to choose from

flat

Selects a volume with a flat top with height capHeight, i.e. a cylinder.

level

Selects a volume whose top surface lies at constant vertical elevation. The value of this elevation is set to capHeight.

Note

In this case only, both capRadius and capHeight must be set.

para (default)

Selects a parabolic volume with formula

\(H(\mathbf{x}) = C_H (1 - |\mathbf{x} - \mathbf{x}_c|^2/C_R)\),

where \(C_H\) is capHeight, \(C_R\) is capRadius and \(\mathbf{x}_c\) is the centre point of the volume’s base at capX, capY.

Cube

A Cube block places a volume of flowing material with a rectangular base into the simulation at the initial time. The following declarations are required:

The location of the volume can be specified by giving either

• the latitude (cubeLat) and longitude (cubeLon) of its centre

or

• the offset of the volume’s centre from the centre of the domain (cubeX, cubeY), in metres.

These required specifiers give:

cubeLat = 13.248745

The latitude of the centre of the volume in decimal degrees in WGS84 coordinates.

cubeLat = 123.686939

The longitude of the centre of the volume in decimal degrees in WGS84 coordinates.

cubeX = 100

The offset of the centre of the volume along the \(x\) axis in metres from the centre of the domain.

cubeY = -50

The offset of the centre of the volume along the \(y\) axis in metres from the centre of the domain.

To set the volume of the release, the following settings are required:

cubeHeight

Sets the height of the volume in metres. N.b. the precise meaning of this depends on the value of cubeShape (see below).

cubeLength

Sets the extent of the volume in the \(x\)-direction.

cubeWidth

Sets the extent of the volume in the \(y\)-direction.

Note

Both cubeLength and cubeWidth should be large enough to ensure that the volume can be represented on the numerical grid.

The following optional parameters set additional attributes for the volume release:

cubeU = 0

Sets the initial velocity \(\bar{u}\) of the release in the \(x\)-direction (m/s).

cubeV = 0

Sets the initial velocity \(\bar{v}\) of the release in the \(y\)-direction (m/s).

cubeConc = 0

Sets the volumetric fraction \(\bar{\psi}\) of solids present in the release.

cubeShape = flat

This selects the geometry of the volume. There are three options to choose from

flat (default)

Selects a volume with a flat top with height cubeHeight, i.e. a cuboid.

level

Selects a volume whose top surface lies at constant vertical elevation. The value of this elevation is set to capHeight.

Source

A Source block specifies conditions for a release of material onto the domain through a time series (referred to as a flux source). A source block is identified using the block keyword Source:.

Multiple flux sources can be added through additional Source blocks.

The flux source is modelled as a circular area through which material is added to the domain at a specified volumetric flux and with a specified solids fraction. The flux source requires a location, size and time series for the volumetric flux and solids fraction.

The location of the source can be specified by giving either

• the latitude (sourceLat) and longitude (sourceLon) of the centre of the source;

or

• the offset of the source centre from the centre of the domain (sourceX, sourceY), in metres.

Note

If using an artificial analytical topographic surface, the location must be set using sourceX, sourceY.

These required specifiers give:

sourceLat = 13.248745

The latitude of the centre of the flux source in decimal degrees in WGS84 coordinates.

sourceLon = 123.686939

The longitude of the centre of the flux source in decimal degrees in WGS84 coordinates.

sourceX = 100

The offset of the centre of the flux source along the \(x\) axis in metres from the centre of the domain.

sourceY = -50

The offset of the centre of the flux source along the \(y\) axis in metres from the centre of the domain.

The following are the additional required settings for a source block:

sourceRadius = 5

The radius of the circular flux source, in metres.

Note

The radius should be large enough to ensure that the source can be represented on the numerical grid.

sourceTime = (  0, 360, 720)

A list of times for which the volumetric flux and solids fraction are given. This takes the form sourceTime = (t0, t1, t2, ..., tN) with ascending times and can contain as many increments as needed.

sourceFlux = (5.0, 7.0, 0.0)

A list of the volume flux (m³/s) at the times given in sourceTime, and takes the form sourceFlux = (Q0, Q1, Q2, ..., QN).

sourceConc = (0.0, 0.0, 0.0)

A list of the volumetric solids concentration values at the times given in sourceTime, and takes the form sourceConc = (psi0, psi1, psi2, ..., psiN).

Note

Each of sourceTime, sourceFlux and sourceConc must contain the same number of points.

Note

Outside the given time series (i.e. for t < t0 or t > tN), Q = 0 and psi = 0.

Between each pair of given time series points, the source flux and concentration values are linearly interpolated.

Output

The Output block sets up output from Kestrel. It is identified using the block keyword Output:.

The only required setting in the Output block is:

N out = 10

The number of output files to be produced. These record simulation data at equally spaced time intervals over the duration of the simulation.

Note

The initial condition is also recorded and does not count towards this total.

The optional settings of the output block are:

base path = ./

A path to a base directory to hold the output directory. Default is current working directory. This is created if it does not exist, providing that permissions allow.

directory = results/

A directory to store the results. This is created if it does not exist and permissions allow.

format = txt

Format of the output files. Options are:

txt (default) – column-headed, comma-delimited text files.

nc or netcdf – NetCDF files. Requires compilation with NetCDF4 (see Installation)

kml – KML files. Requires simulation on georeferenced topography.

Note

Multiple output formats can be specified as a comma-separated list (e.g. format = txt, nc).

Note

KML output feature is currently limited.

info filename = RunInfo.txt

Name of a text file to contain data on the simulation.

The following conditionally optional setting can be given if the format value includes txt, nc or netcdf:

maximums filename = Maximums

The name of a file to contain aggregated data over the duration of the simulation. For a text file, the extension .txt is added. For a NetCDF file, the extension .nc is added.

The following conditionally optional setting can be given if the format value includes txt:

compression = off

Compress text files using tar.

The following conditionally optional setting can be given if format value includes kml:

kml height = 0.1

A threshold depth, in metres, for writing flow depth and flow speed outputs to KML files. Cells with flow depths below kml height are ignored when writing to KML.

inundation time filename = InundationTime

The name of a KML file to store the first time of inundation of points in the domain.

max height filename = MaxHeights

The name of a KML file to store the maximum flow depth, and the time of this maximum, for points in the domain.

max speed filename = MaxSpeeds

The name of a KML file to store the maximum flow speed, and the time of this maximum, for points in the domain.

Parameters

The Parameters block specifies the model closures and associated parameters for the simulation. It is identified using the block keyword Parameters:.

There are two required settings in the Parameters block:

drag = chezy

This sets the basal drag closure function \(\tau_b\), which we assume to be of the form \(\tau_b = \frac{\bar{\rho} \bar{\mathbf{u}}}{\mathbf{u}}\mathcal{F}\), where (as detailed in Physical model) \(\mathcal{F}\) is a friction function that may depend on the local flow variables and any conditionally optional parameters relevant to the selected drag law. Options are:

chezy

Sets the Chézy drag law

\(\mathcal{F} = C_d |\mathbf{u}|^2\)

for modelling the basal resistance of a turbulent fluid.

The conditionally optional setting

chezy co = 0.01

is used to set the value of the constant drag coefficient \(C_d\). It is required to be strictly positive.

coulomb

Sets the Coulomb drag law

\(\mathcal{F} = \mu g_\perp H\)

for modelling granular friction. Here, \(g_\perp\) denotes gravitational acceleration resolved perpendicular to the local slope.

The conditionally optional setting

coulomb co = 0.1

is used to set the value of the constant friction coefficient \(\mu\). It is required to be strictly positive.

manning

Sets the Manning drag law

\(\mathcal{F} = \frac{g_\perp n^2}{H^{1/3}}\),

where \(g_\perp\) is gravitational acceleration resolved perpendicular to the local slope.

The conditionally optional setting

manning co = 0.03

is used to set the dimensional constant coefficient \(n\). It is required to be strictly positive.

pouliquen

Sets Pouliquen and Forterre’s granular drag law (see e.g. Pouliquen & Forterre, J. Fluid Mech., 453, 2002)

\(\mathcal{F} = \mu(I) g_\perp H\),

where \(g_\perp\) is gravitational acceleration resolved perpendicular to the local slope and \(I\) is a dimensionless inertial number, defined by \(I = d |\mathbf{u}| / \sqrt{g_\perp H^3}\). (The parameter \(d\) is the characteristic sediment diameter, which may be set via the solid diameter option.) The friction coefficient is given by

\(\mu(I) = \mu_1 + \frac{\mu_2 - \mu_1}{1 + \beta / I}\)

The following conditionally optional settings are used to set the various empirical coefficients in this model:

pouliquen min = 0.1

Sets \(\mu_1\), which is required to be strictly positive.

pouliquen max = 0.4

Sets \(\mu_2\), which is required to be strictly positive.

pouliquen beta = 0.136

Sets \(\beta\), which is required to be strictly positive.

variable

Sets the following drag law, which interpolates between the Chézy and Pouliquen laws, depending on the solids fraction:

\(\mathcal{F} = (1 - f(\bar{\psi})) C_d |\mathbf{u}|^2 + f(\bar{\psi}) \mu(I) g_\perp H\),

where \(f\) is a switching function, equal to

\(f(\bar{\psi})=\frac{1}{2}[1+\tanh(V_R(\bar{\psi}-\psi^*))]\).

Its parameters may be set via the conditionally optional statements

• voellmy switch rate = 3.0, which sets \(V_R\)

• voellmy switch value = 0.2, which sets \(\psi^*\)

voellmy

Sets the Voellmy drag closure, which is the sum of Chézy and Coulomb drag:

\(\mathcal{F} = C_d |\mathbf{u}|^2 + \mu g_\perp H\).

As above, \(C_d\) and \(\mu\) may be set by the conditionally optional statements chezy co and coulomb co respectively.

erosion = on

This sets the closure function \(\mathcal{E}\) for erosion (see Physical model. If set to anything other than erosion = off, it also activates Kestrel’s morphodynamic capabilities. The following options are available:

fluid

Sets a ‘fluid-like’ erosion, with

\(\mathcal{E} = \max\{ \varepsilon_f u_p (\theta - \theta_c), 0\}\),

where \(u_p = \sqrt{g_\perp d (\rho_s/\rho_f - 1)}\) with \(\rho_s, \rho_f\) denoting sediment and fluid densities respectively. \(\theta = C_d |\mathbf{u}|^2/u_p^2\) denotes the Shields number for the flow. (N.b. \(g_\perp\) and \(d\) are defined in the discussion of the drag setting.)

Erosion occurs when \(\theta\) exceeds the critical value \(\theta_c\), determined by the empirical closure

\(\theta_c = \frac{0.3}{1 + 1.2 R} + 0.055[1-\exp(-0.02R)]\),

where \(R = d [g (\rho_s/\rho_f - 1) / \nu_w^2]^{1/3}\) and \(\nu_w = 1.2\times 10^{-6}\textrm{m}^2/\textrm{s}\) is the kinematic viscosity of water. (cf. Soulsby, Dynamics of Marine Sands, 1997)

Three conditionally optional statements affect this closure:

• The constant coefficient \(\varepsilon_f\) may be set via the

  erosion rate = 1e-3

  whose value is required to be strictly positive.

• Sediment and fluid densities may be set via

  rhos = 2000

  and

  rhow = 1000

  respectively. These are required to be strictly positive.

granular

Sets a ‘granular-like’ erosion, with

\(\mathcal{E} = \max\{ \varepsilon_g u_p [\mu(I)g_\perp H - \mu_n], 0\}\),

where \(\mu(I)\) is Pouliquen’s friction coefficient (see the drag discussion above) and

\(\mu_n = \mu_1 + \frac{\mu_s - \mu_1}{1 + \left(\frac{H}{25d}\right)^2}\)

with \(\mu_s = [\mu_1 + \tan(1^\circ)] / [1 - \mu_1 \tan(1^\circ)]\) (\(\mu_1\) is set by pouliquen min).

The conditionally optional declaration

granular erosion rate = 1e-3

may be used to set the constant coefficient \(\varepsilon_g\), which is required to be strictly positive.

mixed or on (default)

Sets the following erosion law that switches between fluid-like and granular-like erosion rates, depending on the solids fraction:

\(\mathcal{E} = (1 - f(\bar{\psi})) \mathcal{E}_f + f(\bar{\psi}) \mathcal{E}_g\),

where \(\mathcal{E}_f\) and \(\mathcal{E}_g\) are the corresponding erosion rates according to the fluid and granular closures respectively. The function \(f\) is the same switching function as in the case of drag = variable.

off

Sets \(\mathcal{E} = 0\).

simple

Sets a simple model for erosion based on the Shields number, with no critical value:

\(\mathcal{E} = \varepsilon u_p \theta\).

The constant coefficient \(\varepsilon\) may be defined via the conditionally optional setting

erosion rate = 1e-3

which is required to be strictly positive.

The remaining settings in the Parameters block are optional. We list them below:

bed porosity = 0.35

Sets the bed porosity \(p\). This is related to the solid fraction \(\psi_b\) of the bed by \(\psi_b = 1 - p\) and as such, affects the rate of sediment transfer between the flow and bed (see Physical model). Kestrel requires \(0 < p \leq 1\).

Note

In most cases, it is prudent to have bed porosity equal to 1 - maxPack.

deposition = Spearman Manning

This sets the deposition rate closure \(\mathcal{D}\). The following options are available:

none

Sets \(\mathcal{D} = 0\).

simple

Sets a simple quadratic hindered settling law of the form

\(\mathcal{D} = w_s \bar{\psi}(1 - \bar{\psi}/\psi_{\max})\),

where \(w_s\) is characteristic sediment settling speed and \(\psi_{\max}\) is the maximum volume fraction that the flowing sediment may be packed into. These constant coefficients may be set via the conditionally optional declarations:

settling speed = 1e-3

sets \(w_s\). If not explicitly set, Kestrel defaults to using an empirical law based on the solid diameter \(d\):

\(w_s = \frac{\nu_w}{d}\left\{\sqrt{10.36^2 + 1.048R} - 10.36\right\}\),

where (as above) \(R = d [g (\rho_s/\rho_f - 1) / \nu_w^2]^{1/3}\) and \(\nu_w = 1.2\times 10^{-6}\textrm{m}^2/\textrm{s}\) is the kinematic viscosity of water.

maxPack = 0.65

sets \(\psi_{\max}\).

Spearman Manning

Sets an empirical hindered settling law due to Spearman & Manning (Ocean Dynam. 67(3), 2017):

\(\mathcal{D} = w_s \bar{\psi} (1 - \bar{\psi})^a (1 - \bar{\psi}/\psi_{\max})^b\)

The exponents \(a\) and \(b\) are determined via the formulae \(a = 2.7 - 0.15 n\) and \(b = 0.62n - 1.46\) where

\(n = \frac{4.7 + 0.41 (u_p d / \nu_w)^{3/4}}{1 + 0.175 (u_p d / \nu_w)^{3/4}}\)

(Rowe, P. N, Chem. Eng. Sci, 42, 1987). The constants \(w_s\) and \(\psi_{\max}\) may be set via settling speed and maxPack resp. (as above).

eddy viscosity = 0.0

Sets the (constant) value of eddy viscosity \(\nu\) in the model (see Physical model). This value is required to be non-negative.

Warning

If you want to simulate morphodynamics then eddy viscosity must be non-zero. Otherwise, the underlying governing equations are ill-posed as an initial value problem and Kestrel’s numerical solutions will fail to converge as the grid resolution is refined.

erosion critical height = 0.01

Sets a critical flow depth \(H_c\) in metres, below which erosion is not permitted. This ensures that rapid thin flows do not unphysically erode the bed. It is recommended that this is at least equal to the characteristic solid diameter \(d\). It is required to be strictly positive.

A phenomenological function \(\chi(H)\) is pre-multiplied to the value of the morphodynamic transfer rates \(\mathcal{E}\) and \(\mathcal{D}\) to achieve this. It is user-selectable via the conditionally optional setting

morphodynamic damping

This has options

• off. Sets \(\chi = 1\).

• rat3. Sets \(\chi = 0\) if \(H < H_c\), \(\chi = 1\) if \(H > 2 H_c\), \(\chi = (\frac{H}{H_c} - 1)^3 [(2 - \frac{H}{H_c})^{3} + (\frac{H}{H_c} - 1)]^{-1}\) otherwise.

• tanh (default). Sets \(\chi = \frac{1}{2}[1 + \tanh(10 \log(H/H_c))]\).

erosion depth = 1

Sets the depth in metres up to which erosion is permitted. In the notation of Physical model, this means setting \(\Delta b_{\max}\). It is required to be non-negative.

A phenomenological function \(\Theta(\Delta b)\) (where \(\Delta b \equiv b(\mathbf{x},t) - b_0(\mathbf{x},0)\)) is pre-multiplied to the value of the erosion rate \(\mathcal{E}\) to achieve this. It is user-selectable via the conditionally optional setting

erosion transition

This has options

• off. Sets \(\Theta = 1\).

• smooth (default). Sets \(\Theta = \frac{1}{2}[1 + \tanh(10^5(\Delta b + \Delta b_{\max}))]\).

• step. Sets \(\Theta = 0\) if \(\Delta b < -\Delta b_{\max}\), \(\Theta = 1\) otherwise.

g = 9.81

Sets the gravitational acceleration \(g\).

geometric factors = on

This option selects whether the model equations that Kestrel solves should consider geometric corrections that arise due to variations in the topographic surface. This is the default model, described in Physical model. If geometric factors = off, the only effect of slope variation left in the model is that of gravitational forcing along the direction of steepest descent. This is equivalent to setting \(\gamma \equiv 1\) in Physical model.

For more information, see this reference.

Warning

Note that, if geometric factors = off, then \(g_\perp \equiv g\).

solid diameter = 1e-3

Sets the characteristic sediment diameter \(d\). This affects various optional closures, that model the physics of grains such as the settling speed, Pouliquen drag rule and granular erosion. It is required to be strictly positive.

Solver

The Solver block sets up the core numerical solver in Kestrel. The solver block is identified using the block keyword Solver:.

The only required setting in the solver block is:

T end

The time, in seconds, at which to end the simulation.

The optional settings of the solver block are:

limiter = MinMod2

The slope limiter using to approximate gradients in the governing equations. Options are:

MinMod1

The minmod limiter, defined as

\[\mathrm{MinMod}(a, b) = \tfrac{1}{2}(\mathrm{sgn}(a) + \mathrm{sgn}(b)).\min(\left|a\right|, \left|b\right|).\]

for real numbers \(a\) and \(b\) and \(\mathrm{sgn}(a)\) denotes the signum function.

MinMod2 (default)

The generalised minmod limiter, defined as

\[\begin{split}\mathrm{MinMod}(a, b) = \begin{cases} \mathrm{sgn}(a) . \min(\theta\left|a\right|, \tfrac{1}{2}\left|a+b\right|, \theta\left|b\right|), & \text{if $\mathrm{sgn}(a)=\mathrm{sgn}(b)$}\\ 0, & \text{if $\mathrm{sgn}(a)+\mathrm{sgn}(b)$ = 0}. \end{cases}\end{split}\]

The parameter \(1\le \theta \le 2\) and we impose \(\theta = 1.3\) which we have found to work well (see also Chertock et al. (2015)).

Note

MinMod2 is strongly recommended.

van albada or albada

The van Albada limiter is a smooth function that tends to a central difference approximation in smooth regions, and across discontinuities the slope is biased to smallest value of the two one-sided slopes. In general, the van Albada limiter has a parameter, but here we use a version with parameter set to zero (see Lu & Tadmor (2023))

\[\mathrm{vanAlbada}(a,b) = \frac{a^2 b + a b^2}{a^2 + b^2}.\]

WENO

A weighted essentially non-oscillatory type limiter.

\[\mathrm{WENO}(a,b) = \frac{w(a)a + w(b)b}{w(a) + w(b)} \text{ with weighting } w(x) = (x^2+\epsilon)^{-2}\]

Kestrel imposes \(\epsilon = 10^{-6}\).

None

No slope limiter used. Gradients are approximated using central differences.

Note

As the governing equations are hyperbolic, solutions can exhibit discontinuities. Therefore, None is not recommended.

height threshold = 1e-6

A threshold on flow depths. Flow depths below height threshold are neglected.

Tile buffer = 1

Neighbouring tiles are activated if flow reaches tile buffer cells from a tile edge. The default value tile buffer = 1 should ensure that a neighbouring tile is added when needed.

CFL = 0.5 for 1D simulations; CFL = 0.25 for 2D simulations

The Courant-Friedrichs-Lewy (CFL) number of the simulation. This determines the maximum time step.

• For 1D simulations the scheme requires CFL \(\le 0.5\).

• For 2D simulations the scheme requires CFL \(\le 0.25\).

max dt = HUGE(1.0)

The maximum time step. Note that the default value ensures that the maximum time step is the either the time step determined by the CFL condition, or the output time step.

T start = 0

The start time of the simulation, in seconds.

Restart = off

Should the simulation restart from a previous result? Options are:

on – restart from a previous simulation. This requires an existing Kestrel results directory containing the RunInfo.txt file, which is used to determine suitable simulation parameters and to locate the last output file to be used as the initial condition.

off – start a new simulation.

Note

This feature is useful if a simulation is interrupted for some reason. By selecting on and re-running, the simulation will pick up from where it left off.

Initial condition

Specifies the path to a Kestrel result file to be used as an initial condition. On start-up Kestrel loads the solution fields from this file and simulates forward from this point.

If Restart = on, then RunInfo.txt is used to determine the simulation parameters. Otherwise, (by default) they are read in from the usual input file given on the command line.

The following conditionally optional variables are used only if Boundary conditions = sponge in the Domain block:

Sponge strength = 0.2

When using a sponge layer boundary condition, the solution’s quantities are gradually damped on the tiles bordering the domain boundary. The damping rate is set by the Sponge strength settings.

Note

Care must be taken in setting Sponge strength. If the damping is too weak, flow quantities may be non-zero at the domain edge, causing errors. If the damping is too strong, flow quantities in the interior can be influenced by those in the sponge layer tiles. The flow in the sponge layer tiles should be discarded.

Topog

The Topog block specifies the topography to be used for the simulation. Topography can be given either as a Digital Elevation Model (DEM) or through a parameterised analytic function. The Topog block is identified using the block keyword Topog:.

There is a single required parameter in the Topog block:

Type

The type of topography to use. Options are:

DEM or Raster

Use a georeferenced DEM contained in a geotiff file.

SRTM

Use a georeferenced SRTM file, contained in an SRTM archive.

Function

Use a parameterised analytic function defining a surface.

The following statement is conditionally required when Type = DEM or Type = Raster:

raster file

The name of the georeferenced raster file containing the DEM.

The following conditionally optional parameters may be used when Type = DEM or Type = Raster:

dem directory

The directory containing raster file. Default is the current working directory.

embed raster

Should the raster DEM contained in raster file be embedded in SRTM topography? Options are:

• Yes or On – embed the DEM within SRTM data.

  Warning

  If embed raster = on then flows will transition between the DEM topography to SRTM at the edge of the DEM. While this can be useful for simulations where the available DEM coverage is partial, any offsets between the DEM and SRTM surface maps can cause problems.

• No or Off – do no embed the DEM.

  Warning

  If embed raster = off then error will occur if the flow reaches the edge of the DEM.

The following parameter is conditionally required when Type = DEM or Type = Raster and embed raster = on. Additionally, it is conditionally required when Type = SRTM.

srtm directory

The directory containing the SRTM data. This should be in the form of an SRTM archive, i.e. either: - containing zip files with extension .SRTMGL1.hgt.zip - a directory with a structure N01/N01*.tif, N02/N02*.tif, …

Note

SRTM files may be either geotiff or USGS hgt format (including zipped hgt files with extension .SRTMGL1.hgt.zip).

There is one conditionally required setting needed when Type = Function – Topog function giving the choice of function \(b_0\) to be used for the initial bed elevation. This is accompanied by a conditionally optional setting – Topog params, which specifies a list of parameters to use with the function. The Topog params variable takes the form Topog params = (a, b, c, ...). Its values convey a particular meaning for each function choice, as described below for the currently implemented functions:

Topog function = flat

A flat surface, \(b_0(x,y) = 0\).

Topog params not used.

Topog function = xslope

A constant slope along the \(x\)-coordinate, \(b_0(x,y) = \alpha x\).

Topog params = (alpha)

• alpha – The slope of the plane.

Topog function = yslope

A constant slope along the \(y\)-coordinate, \(b_0(x,y) = \beta y\).

Topog params = (beta)

• beta – The slope of the plane.

Topog function = xyslope

A general plane inclined along both the \(x\)- and \(y\)-coordinates, \(b_0(x,y) = \alpha x + \beta y\).

Topog params = (alpha, beta)

• alpha – The slope of the plane along \(x\).

• beta – The slope of the plane along \(y\).

Topog function = x2slopes

A surface with two slopes in the \(x\)-direction, connected by a circular arc near \(x = 0\).

Topog params = (alpha, beta, R)

• alpha – slope on left-hand-side.

• beta – slope on right-hand-side.

• R – radius of connecting circular arc.

Topog function = xBislope

A surface with two slopes in the limits \(x\to\pm\infty\), connected by a smooth transition. The surface has the form

\(b_{0}(x,y) = -\tfrac{1}{2}\left(\tan\phi_{1} + \tan\phi_{2}\right)x + \tfrac{1}{2}\left(\tan\phi_{1} - \tan\phi_{2}\right)\lambda\log\left[\cosh\left(x/\lambda\right)\right].\)

Topog params = (phi1, phi2, lambda)

• phi1 – the slope angle for \(x\to -\infty\), in degrees. A positive value corresponds to an elevation decreasing from left to right.

• phi2 – the slope angle for \(x\to +\infty\), in degrees. A positive value corresponds to an elevation decreasing from left to right.

• lambda – the characteristic length scale of the smooth transition region.

Topog function = USGS

Parametrisation of the USGS flume. This has slope of 31° for \(x<0\), and slope 2.4° for \(x>x_{1}>0\) that are connected by a smooth \(\cosh\) curve section. Note \(x_{1}\) is determined to ensure smooth connection of the slope elements. The flume is confined by walls for \(x<8.5\) m, that are represented as \(\tanh\) profile humps. See Iverson et al. (2010) for details.

Topog params = (wallH, sigma)

• wallH – the height of the sidewalls of the flume.

• sigma – the width of the sidewalls of the flume.

Topog function = xsinslope

One-dimensional sinusoidal variation along the x-direction, with one complete period in the specified domain. Letting \(L_{x}\) denote the domain length in \(x\), the surface is \(b_{0}(x,y) = \epsilon \sin(2\pi x / L_{x}).\)

Topog params = (epsilon)

• epsilon – the amplitude of the sinusoidal variation.

Topog function = xysinslope

Two-dimensional sinusoidal variation, with one complete period in the specified domain. Letting \(L_{x}\) and \(L_{y}\) denote the domain lengths in \(x\) and \(y\) respectively, the surface is \(b_{0}(x,y) = \epsilon \sin(2\pi x / L_{x}) \sin(2\pi y / L_{y}).\)

Topog params = (epsilon)

• epsilon – the amplitude of the sinusoidal variation.

Topog function = xhump

One-dimensional cosine hump on a flat topography,

\(b_{0}(x,y) = \tfrac{1}{2} A \left(1 + \cos(\pi x/L)\right)\),

for \(-L \le x \le L\).

Topog params = (A, L)

• A – the amplitude of the hump.

• L – the half-length of the hump.

Topog function = xtanh

One-dimensional \(\tanh\) surface,

\(b_{0}(x,y) = A\left[ 1 + \tanh\left((x-x_{0})/L\right) \right]\)

Topog params = (x0, A, L)

• x0 – the centre of the tanh profile.

• A – the amplitude of the hump.

• L – the half-length of the hump.

Topog function = xparab

One-dimensional parabolic surface,

\(b_{0}(x,y) = Ax^{2}\)

Topog params = (A)

• A – coefficient of the parabola.

Topog function = xyparab

Two-dimensional parabolic surface,

\(b_{0}(x,y) = Ax^{2} + By^{2}\)

Topog params = (A, B)

• A – coefficient of \(x^{2}\) for the parabola.

• B – coefficient of \(y^{2}\) for the parabola.