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
andnYpertile = 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 variablesBoundary Hn
,Boundary U
,Boundary V
, andBoundary 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
andcapHeight
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\) iscapRadius
and \(\mathbf{x}_c\) is the centre point of the volume’s base atcapX
,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
andcubeWidth
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 (m3/s) at the times given in
sourceTime
, and takes the formsourceFlux = (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 formsourceConc = (psi0, psi1, psi2, ..., psiN)
.Note
Each of
sourceTime
,sourceFlux
andsourceConc
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
ornetcdf
– 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
andcoulomb 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
oron
(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
andgranular
closures respectively. The function \(f\) is the same switching function as in the case ofdrag = 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 to1 - 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
andmaxPack
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 andgranular
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
oralbada
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 valuetile buffer = 1
should ensure that a neighbouring tile is added when needed.
CFL = 0.5
for 1D simulations;CFL = 0.25
for 2D simulationsThe 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
orRaster
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
orOn
– 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
orOff
– 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.