Model in- & output
Setting up an AeoLiS model involves configuring various parameters and input files. In this section, we will cover:
Model Input: The input files and their formats.
Model Output: How to analyse the output.
Default settings: Overview of the model parameters, incl. default values.
Activate/deactivate processes: Process flags.
Model state/output: Overview of model state variables.
Guidance on solver use
A general tip for setting up an AeoLiS model is to start simple and build up the complexity. As a starting point, use the default values and then slowly deviate from those and turn on processes. Start with a relatively coarse grid to allow for faster simulation testing. It is also highly recommended to start with constant conditions, such as a one-directional constant wind, making it easier to interpret if the model output is logical.
While testing, go through the sequential steps AeoLiS takes, plot the variables, and critically assess whether the results are expected. Important model state variables to check include: wind speed and direction (uws, uwn), shear velocity (ustars, ustarn), velocity threshold (uth), sediment concentrations (Cu, Ct), pickup (pickup), and, finally, the resulting bed level change (zb, dzb).
In case you run into issues, we encourage users to post questions and case studies on the AeoLiS Discussion Board. We use this public forum so our help and advice is available to everyone.
Model input
The computational grid and boundary conditions for AeoLiS are specified through external
input files called by the main configuration file. The computational domain is defined
using specific spatial files (*.grd), while boundary conditions for wind, wave, and tides
are provided via *.txt files. An overview of these files is provided in the table below.
Input File |
Keyword |
Dimensions |
Requirement |
File Description |
|---|---|---|---|---|
aeolis.txt |
N/A |
N/A |
Mandatory |
Main file containing parameter definitions |
x.grd |
|
(ny, nx) |
Mandatory |
File containing cross-shore grid coordinates |
y.grd |
|
(ny, nx) |
Mandatory |
File containing alongshore grid coordinates |
z.grd |
|
(ny, nx) |
Mandatory |
File containing topography and bathymetry data (bed level) |
zne.grd |
|
(ny, nx) |
Optional |
File containing the non-erodible layer elevation |
veg.grd |
|
(ny, nx) |
Optional |
Initial vegetation density (if |
hveg.grd |
|
(ny*nx*nspecies) |
Optional |
Vegetation height per species (if |
Nt.grd |
|
(ny*nx*nspecies) |
Optional |
Vegetation density/tillers per species (if |
mass.txt |
|
(nx*ny, nfractions*nlayers) |
Optional |
Sediment mass data (for space-varying grain sizes) |
wind.txt |
|
(ntimesteps, 3) |
Mandatory |
File containing wind speed and direction data |
tide.txt |
|
(ntimesteps, 2) |
Optional |
Water elevation data (if |
wave.txt |
|
(ntimesteps, 3) |
Optional |
Wave height and period data (if |
meteo.txt |
|
(ntimesteps, 6) |
Optional |
Meteorological data (if |
mask_*.grd |
|
(ny, nx) |
Optional |
Spatial masks for tide, wave, or runup boundary conditions |
Configuration file (aeolis.txt)
The main configuration file controls all processes and input files.
Parameters in the file are specified by various keywords; each keyword has a pre-defined
default value (see constants.py and the Default settings tab) that will be used if it is not directly specified. Among the keywords
are those defining grid files (xgrid_file, ygrid_file,
bed_file) and boundary conditions (tide_file, wave_file, wind_file).
Physical processes in AeoLiS can be toggled by setting process keywords to True (T) or False (F).
Example parameter files can be found in the examples folder on the AeoLiS GitHub.
Grid files (*.grd)
All grid files (x.grd, y.grd, z.grd, zne.grd, veg.grd, etc.) must have the exact same dimensions.
Each value (or element) within these matrices represents a single computational cell. This means that element [i, j]
in the x.grd file corresponds to the exact same physical cell as element [i, j] in the z.grd file.
The model can be run in 1D and 2D mode depending on the input dimensions. Most processes are implemented in 2D mode. To run the model in 2D mode, all grid files should contain 2D matrices (ny, nx) of the same size. To run the model in 1D mode, all grid files should contain 1D vectors (nx, 1). Because some processes are easier to solve in 2D, the 1D model is internally converted to a quasi-2D model by repeating the vectors three times (nx, 3), assuming the same resolution spacing as the cross-shore direction. Results are then converted back to 1D by extracting the middle vector.
x.grd and y.grd (coordinates)
The x.grd and y.grd files define the coordinates of the computational grid in meters. Important: The structure and order of the elements in these files determines the boundary definitions (i.e., which boundary is on which edge, what is longshore, what is cross-shore), while the specific coordinates determine the location and orientation of the domain. To ensure the model correctly interprets model grids, generation should follow this step-by-step approach (see also Python example further down below):
Define dimensions and resolution: Set your domain length and grid spacing. At this stage, the x-direction always represents cross-shore and y-direction longshore. The model does not (yet) support variable grid sizes; the cross-shore resolution must equal the alongshore resolution (e.g.,
dx = 10; dy = 10).Create domain axes: Generate 1D arrays for the cross-shore (x) and alongshore (y) directions. Both arrays must be ascending. In case of the x-array, the first element (
x[0]) is at the onshore boundary, and the last element (x[-1]) is the offshore boundary (e.g.,x = np.arange(0, Lx + dx, dx)).Generate the grid: Use the 1D axes to create 2D coordinate matrices (e.g.,
X, Y = np.meshgrid(x, y)).Shift and rotate: Now, shift and rotate the coordinatet to match the desired real-world location and orientation (e.g.,
Xr = X * np.cos(theta) - Y * np.sin(theta) + x0). Following this order ensures both the structure of the files and the actual coordinates are correct.
By handling the physical orientation through the grid coordinates, the model orientation keyword alfa is no longer required and should be left untouched.
z.grd (bed level)
The z.grd file provides the surface elevation for every cell defined in the spatial grids. Elevation values should be defined such that positive values are above the vertical datum (up) and negative values are below (down).
zne.grd (non-erodible Layer)
The zne.grd file defines the elevation of the non-erodible layer beneath the sand surface. The model will not erode sediment below this elevation. It follows the same dimensions and elevation datum as z.grd.
veg.grd (vegetation)
The veg.grd file is an optional grid providing the initial vegetation coverage (density) at each cell (0-1).
Fig. 12 File format for a 1D AeoLis vegetation grid. Each red dot is the vegetation density at a specific location in the computational grid.
Masks (tide, wave, and runup)
Masks (e.g., mask_tide, mask_wave, mask_runup) can be used to spatially modify boundary conditions. This can be useful when, for instance, onshore elevations are lower than the offshore water level but remain dry, or when water bodies are disconnected from the offshore (such as a barrier island where the lagoon has no tidal or wave action).
Masks use complex numbers to apply both a scaling multiplier (real) and a static offset (complex) to the boundary condition. For instance, for waveheight (Hs): Hs = $\mathbb{R}$ * Hs + $\mathbb{C}$. To half the boundary value in a specific area (e.g., sheltered waves inland), use 0.5. To set a fixed static value regardless of the incoming boundary condition (e.g., an inland lake permanently fixed at +2m water level), use a complex value such as 0 + 2j.
Example Python script
Generating the grid files described above
import numpy as np
# 1. Define dimensions and resolution
Lx, Ly = 1000, 500 # Lx is always cross-shore, Ly is longshore
dx = 10 # dx must equal dy
x0, y0 = 5000, 5000 # coordinates for a grid corner point
theta = np.radians(30) # grid rotation in radians
# 2. Create axes
x = np.arange(0, Lx + dx, dx) # must be ascending (first=onshore, last=offshore)
y = np.arange(0, Ly + dx, dx)
# 3. Create orthogonal grid
X, Y = np.meshgrid(x, y)
# 4. Rotate and shift domain to real-world coordinates
Xr = X * np.cos(theta) - Y * np.sin(theta) + x0
Yr = X * np.sin(theta) + Y * np.cos(theta) + y0
# *Because orientation is baked into the coordinates, set `alfa = 0` in aeolis.txt*
# 5. Populate specific grid variables
z = np.zeros_like(X) # example: Flat bed level at 0.0m
zne = np.full_like(z, -2.0) # example: Flat non-erodible layer at -2.0m
veg = np.zeros_like(z) # example: no initial vegetation
# 6. Create masks (e.g., tide mask for an inland lake)
tide_mask = np.ones_like(X, dtype=complex)
tide_mask[:, :10] = 0.0 + 2.0j # 0.0 multiplier (no tidal variation) + 2.0j elevated
# Save outputs to text files (to be read as .grd)
np.savetxt('x.grd', Xr)
np.savetxt('y.grd', Yr)
np.savetxt('z.grd', z)
np.savetxt('zne.grd', zne)
np.savetxt('veg.grd', veg)
np.savetxt('mask_tide.grd', tide_mask)
Multi-dimensional input
Some variables require multiple dimensions per spatial grid cell, such as spatially varying grain sizes (multiple fractions and bed layers) or complex vegetation cover (multiple species). AeoLiS requires these multi-dimensional arrays to be flattened into lower-dimensional formats before they can be saved as textfiles.
mass.txt: Defines the mass of each sediment fraction per bed layer. If the grain size distribution is uniform across the domain, which is most often the case, this file is not needed (see keywords
grain_distandgrain_sizein the configuration file). The 4D shape of (ny, nx, nlayers, nfractions) must be reshaped into a 2D matrix of shape (nx * ny, nfractions * nlayers). The rows represent the flattened spatial coordinates, and the columns are grouped by bed layer (e.g., Layer 1: Fraction 1, Fraction 2; Layer 2: Fraction 1, Fraction 2).hveg.grd and Nt.grd: When using the grass vegetation method (
method_vegetation = grass), the model needs vegetation height (hveg) and density (Nt). Because this method supports multiple interacting species, the 3D shape of (ny, nx, nspecies) is flattened into a 1D vector of shape (ny * nx * nspecies).
Fig. 13 File format for a 2D AeoLis mass file for spatially variable grain size distributions. Each red dot is the mass for a specific sediment fraction within a specific bed layer at a given spatial location.
Example Python script
Loading spatial grid dimensions and generating multi-dimensional input grids.
import numpy as np
# 1. Load spatial grid to get dimensions
X = np.loadtxt('x.grd')
ny, nx = X.shape
n_cells = ny * nx
n_species = 2
# 2a. Create mass.txt (sediment fractions and layers)
n_layers = 3
n_fractions = 4
mass_per_fraction = np.array([0.4, 0.3, 0.2, 0.1]) # kg per fraction (uniform sediment distribution across the domain)
# Create the mass matrix using tile
mass_matrix = np.tile(mass_per_fraction, (n_cells, n_layers)) # (n_cells, n_fractions * n_layers)
# Initialize the vegetation arrays
hveg = np.zeros((ny, nx, n_species))
Nt = np.zeros((ny, nx, n_species))
# Fill the vegetation arrays
hveg[:, :, 0] = 0.5 # first species is 0.5m tall
hveg[:, :, 1] = 1.0 # second species is 1.0m tall
Nt[:, :, 0] = 10.0 # tillers/m2
Nt[:, :, 1] = 5.0
# Flatten the vegetation arrays
hveg_flat = hveg.reshape(-1)
Nt_flat = Nt.reshape(-1)
# Save files
np.savetxt('mass.txt', mass_matrix)
np.savetxt('hveg.grd', hveg_flat)
np.savetxt('Nt.grd', Nt_flat)
Time-series (*.txt)
Environmental forcing in AeoLiS (wind, water, waves) is provided through time-series files. For all of these files, the format follows the same structure: the first column represents time in seconds w.r.t. refdate in the configruation file. The other columns contain the variables at those given times. The model will automatically interpolate these data points to match the modelling time steps.
wind.txt
The wind.txt file provides the wind boundary conditions driving aeolian transport. It contains three columns: 1. Time (s), 2. Wind speed (m/s), 3. Wind direction (degrees). Wind directions can be specified in either nautical or cartesian convention, which is set using the wind_convention keyword in the configuration file.
Fig. 14 File format for wind boundary conditions file for AeoLis input.
tide.txt
The tide.txt file contains the water elevation data for the duration of the simulation. It contains two columns: 1. Time (s), 2. Water elevation (m)
Fig. 15 File format for the water elevation conditions file for AeoLis input.
wave.txt
The wave.txt file provides the wave data used by AeoLiS to calculate runup. It contains three columns: 1. Time (s), 2. Significant wave height (m), 3. Peak wave period (s)
Fig. 16 File format for the wave conditions file for AeoLis input.
meteo.txt
The meteo.txt file contains meteorological data and is only required if using the groundwater module by Hallin (2023) to simulate surface moisture (process_groundwater = T). It contains six columns: 1. Time (s), 2. Temperature (°C), 3. Precipitation (mm/hr), 4. Relative humidity (%), 5. Global radiation (MJ/$m^2$/day), 6. Air pressure (kPa)
Fig. 17 File format for meteorological data used to simulate surface moisture in AeoLiS.
Example Python script
Generating and saving time-series files. In this example, we generate a 10-day simulation with a constant onshore wind, a harmonic tide, and constant wave conditions.
import numpy as np
# 1. Define time array
days = 10
dt = 3600 # 1-hour intervals
time_sec = np.arange(0, days * 24 * 3600 + dt, dt)
# 2. Wind (constant onshore wind)
# Assuming 270 degrees is straight onshore
wind_speed = np.full_like(time_sec, 8.0) # 8 m/s
wind_dir = np.full_like(time_sec, 270.0) # 270 degrees
wind_data = np.column_stack((time_sec, wind_speed, wind_dir))
# 3. Tide (harmonic tide)
tide_amp = 1.0 # 1 meter amplitude
tide_period = 12 * 3600 # 12-hour period in seconds
water_level = tide_amp * np.sin(2 * np.pi * time_sec / tide_period)
tide_data = np.column_stack((time_sec, water_level))
# 4. Waves (constant wave height and period)
wave_height = np.full_like(time_sec, 1.5) # 1.5m significant wave height
wave_period = np.full_like(time_sec, 6.0) # 6s peak period
wave_data = np.column_stack((time_sec, wave_height, wave_period))
# 5. Save to text files
np.savetxt('wind.txt', wind_data, fmt='%.2f')
np.savetxt('tide.txt', tide_data, fmt='%.2f')
np.savetxt('wave.txt', wave_data, fmt='%.2f')
Visualizating input settings
When running the model, you can automatically generate diagnostic plots by setting the keyword visualization = T in the configuration file. This will automatically generate figures in your simulation folder to help verify your setup:
figure_grid_initialization.png: Displays the grid orientation and boundary definitions.figure_params_initialization.png: Shows the most relevant spatial parameters mapped onto the domain.figure_timeseries_initialization.png: Visualizes the time series of your environmental boundary conditions.
All parameters and defaults
All available configuration parameters, and their default values, are collected in constants.py.
For any configuration parameters not defined in the model configuration file,
the default value is used.
DEFAULT_CONFIG = {
'process_wind' : True, # Enable the process of wind
'process_transport' : True, # Enable the process of transport
'process_bedupdate' : True, # Enable the process of bed updating
'process_threshold' : True, # Enable the process of threshold
'th_grainsize' : True, # Enable wind velocity threshold based on grainsize
'th_bedslope' : False, # Enable wind velocity threshold based on bedslope
'th_moisture' : False, # Enable wind velocity threshold based on moisture
'th_drylayer' : False, # Enable threshold based on drying of layer
'th_humidity' : False, # Enable wind velocity threshold based on humidity
'th_salt' : False, # Enable wind velocity threshold based on salt
'th_sheltering' : False, # Enable wind velocity threshold based on sheltering by roughness elements
'th_nelayer' : False, # Enable wind velocity threshold based on a non-erodible layer
'process_avalanche' : False, # Enable the process of avalanching
'process_shear' : False, # Enable the process of wind shear
'process_tide' : False, # Enable the process of tides
'process_wave' : False, # Enable the process of waves
'process_runup' : False, # Enable the process of wave runup
'process_moist' : False, # Enable the process of moist
'process_mixtoplayer' : False, # Enable the process of mixing
'process_wet_bed_reset' : False, # Enable the process of bed-reset in the intertidal zone
'process_meteo' : False, # Enable the process of meteo
'process_salt' : False, # Enable the process of salt
'process_humidity' : False, # Enable the process of humidity
'process_groundwater' : False, # Enable the process of groundwater
'process_scanning' : False, # Enable the process of scanning curves
'process_inertia' : False, # NEW
'process_separation' : False, # Enable the including of separation bubble
'process_vegetation' : False, # Enable the process of vegetation
'process_fences' : False, # Enable the process of sand fencing
'process_dune_erosion' : False, # Enable the process of wave-driven dune erosion
'process_seepage_face' : False, # Enable the process of groundwater seepage (NB. only applicable to positive beach slopes)
'visualization' : False, # Boolean for visualization of model interpretation before and just after initialization
'output_sedtrails' : False, # NEW! [T/F] Boolean to see whether additional output for SedTRAILS should be generated
'nfraction_sedtrails' : 0, # [-] Index of selected fraction for SedTRAILS (0 if only one fraction)
'xgrid_file' : None, # Filename of ASCII file with x-coordinates of grid cells
'ygrid_file' : None, # Filename of ASCII file with y-coordinates of grid cells
'bed_file' : None, # Filename of ASCII file with bed level heights of grid cells
'wind_file' : None, # Filename of ASCII file with time series of wind velocity and direction
'tide_file' : None, # Filename of ASCII file with time series of water levels
'wave_file' : None, # Filename of ASCII file with time series of wave heights
'meteo_file' : None, # Filename of ASCII file with time series of meteorlogical conditions
'bedcomp_file' : None, # Filename of ASCII file with initial bed composition
'threshold_file' : None, # Filename of ASCII file with shear velocity threshold
'fence_file' : None, # Filename of ASCII file with sand fence location/height (above the bed)
'ne_file' : None, # Filename of ASCII file with non-erodible layer
'veg_file' : None, # Filename of ASCII file with initial vegetation density
'supply_file' : None, # Filename of ASCII file with a manual definition of sediment supply (mainly used in academic cases)
'wave_mask' : None, # Filename of ASCII file with mask for wave height
'tide_mask' : None, # Filename of ASCII file with mask for tidal elevation
'runup_mask' : None, # Filename of ASCII file with mask for run-up
'threshold_mask' : None, # Filename of ASCII file with mask for the shear velocity threshold
'gw_mask' : None, # Filename of ASCII file with mask for the groundwater level
'vver_mask' : None, #NEWBvW # Filename of ASCII file with mask for the vertical vegetation growth
'nx' : 0, # [-] Number of grid cells in x-dimension
'ny' : 0, # [-] Number of grid cells in y-dimension
'dt' : 60., # [s] Time step size
'dx' : 1.,
'dy' : 1.,
'CFL' : 1., # [-] CFL number to determine time step in explicit scheme
'accfac' : 1., # [-] Numerical acceleration factor
'max_bedlevel_change' : 999., # [m] Maximum bedlevel change after one timestep. Next timestep dt will be modified (use 999. if not used)
'tstart' : 0., # [s] Start time of simulation
'tstop' : 3600., # [s] End time of simulation
'restart' : None, # [s] Interval for which to write restart files
'dzb_interval' : 86400, # [s] Interval used for calcuation of vegetation growth
'output_times' : 60., # [s] Output interval in seconds of simulation time
'output_file' : None, # Filename of netCDF4 output file
'output_vars' : ['zb', 'zs',
'Ct', 'Cu',
'uw', 'udir',
'uth', 'mass'
'pickup', 'w'], # Names of spatial grids to be included in output
'output_types' : [], # Names of statistical parameters to be included in output (avg, sum, var, min or max)
'external_vars' : [], # Names of variables that are overwritten by an external (coupling) model, i.e. CoCoNuT
'grain_size' : [225e-6], # [m] Average grain size of each sediment fraction
'grain_dist' : [1.], # [-] Initial distribution of sediment fractions
'nlayers' : 3, # [-] Number of bed layers
'layer_thickness' : .01, # [m] Thickness of bed layers
'g' : 9.81, # [m/s^2] Gravitational constant
'v' : 0.000015, # [m^2/s] Air viscosity
'rhoa' : 1.225, # [kg/m^3] Air density
'rhog' : 2650., # [kg/m^3] Grain density
'rhow' : 1025., # [kg/m^3] Water density
'porosity' : .4, # [-] Sediment porosity
'Aa' : .085, # [-] Constant in formulation for wind velocity threshold based on grain size
'z' : 10., # [m] Measurement height of wind velocity
'h' : None, # [m] Representative height of saltation layer
'k' : 0.001, # [m] Bed roughness
'L' : 100., # [m] Typical length scale of dune feature (perturbation)
'l' : 10., # [m] Inner layer height (perturbation)
'c_b' : 0.2, # [-] Slope at the leeside of the separation bubble # c = 0.2 according to Durán 2010 (Sauermann 2001: c = 0.25 for 14 degrees)
'mu_b' : 30, # [deg] Minimum required slope for the start of flow separation
'buffer_width' : 10, # [m] Width of the bufferzone around the rotational grid for wind perturbation
'sep_filter_iterations' : 0, # [-] Number of filtering iterations on the sep-bubble (0 = no filtering)
'zsep_y_filter' : False, # [-] Boolean for turning on/off the filtering of the separation bubble in y-direction
'Cb' : 1.5, # [-] Constant in bagnold formulation for equilibrium sediment concentration
'Ck' : 2.78, # [-] Constant in kawamura formulation for equilibrium sediment concentration
'Cl' : 6.7, # [-] Constant in lettau formulation for equilibrium sediment concentration
'Cdk' : 5., # [-] Constant in DK formulation for equilibrium sediment concentration
# 'm' : 0.5, # [-] Factor to account for difference between average and maximum shear stress
# 'alpha' : 0.4, # [-] Relation of vertical component of ejection velocity and horizontal velocity difference between impact and ejection
'kappa' : 0.41, # [-] Von Kármán constant
'sigma' : 4.2, # [-] Ratio between basal area and frontal area of roughness elements
'beta' : 130., # [-] Ratio between drag coefficient of roughness elements and bare surface
'bi' : 1., # [-] Bed interaction factor
'T' : 1., # [s] Adaptation time scale in advection equation
'Tdry' : 3600.*1.5, # [s] Adaptation time scale for soil drying
'Tsalt' : 3600.*24.*30., # [s] Adaptation time scale for salinitation
'Tbedreset' : 86400., # [s]
'eps' : 1e-3, # [m] Minimum water depth to consider a cell "flooded"
'gamma' : .5, # [-] Maximum wave height over depth ratio
'xi' : .3, # [-] Surf similarity parameter
'facDOD' : .1, # [-] Ratio between depth of disturbance and local wave height
'csalt' : 35e-3, # [-] Maximum salt concentration in bed surface layer
'cpair' : 1.0035e-3, # [MJ/kg/oC] Specific heat capacity air
'fc' : 0.11, # [-] Moisture content at field capacity (volumetric)
'w1_5' : 0.02, # [-] Moisture content at wilting point (gravimetric)
'resw_moist' : 0.01, # [-] Residual soil moisture content (volumetric)
'satw_moist' : 0.35, # [-] Satiated soil moisture content (volumetric)
'resd_moist' : 0.01, # [-] Residual soil moisture content (volumetric)
'satd_moist' : 0.5, # [-] Satiated soil moisture content (volumetric)
'nw_moist' : 2.3, # [-] Pore-size distribution index in the soil water retention function
'nd_moist' : 4.5, # [-] Pore-size distribution index in the soil water retention function
'mw_moist' : 0.57, # [-] m, van Genucthen param (can be approximated as 1-1/n)
'md_moist' : 0.42, # [-] m, van Genucthen param (can be approximated as 1-1/n)
'alfaw_moist' : -0.070, # [cm^-1] Inverse of the air-entry value for a wetting branch of the soil water retention function (Schmutz, 2014)
'alfad_moist' : -0.035, # [cm^-1] Inverse of the air-entry value for a drying branch of the soil water retention function (Schmutz, 2014)
'thick_moist' : 0.002, # [m] Thickness of surface moisture soil layer
'K_gw' : 0.00078, # [m/s] Hydraulic conductivity (Schmutz, 2014)
'ne_gw' : 0.3, # [-] Effective porosity
'D_gw' : 12, # [m] Aquifer depth
'tfac_gw' : 10, # [-] Reduction factor for time step in ground water calculations
'Cl_gw' : 0.7, # [m] Groundwater overheight due to runup
'in_gw' : 0, # [m] Initial groundwater level
'GW_stat' : 1, # [m] Landward static groundwater boundary (if static boundary is defined)
'max_moist' : 10., # NEWCH # [%] Moisture content (volumetric in percent) above which the threshold shear velocity is set to infinity (no transport, default value Delgado-Fernandez, 2010)
'max_moist' : 10., # [%] Moisture content (volumetric in percent) above which the threshold shear velocity is set to infinity (no transport, default value Delgado-Fernandez, 2010)
'theta_dyn' : 33., # [degrees] Initial Dynamic angle of repose, critical dynamic slope for avalanching
'theta_stat' : 34., # [degrees] Initial Static angle of repose, critical static slope for avalanching
'avg_time' : 86400., # [s] Indication of the time period over which the bed level change is averaged for vegetation growth
'gamma_vegshear' : 16., # [-] Roughness factor for the shear stress reduction by vegetation
'hveg_max' : 1., # [m] Max height of vegetation
'dzb_opt' : 0., # [m/year] Sediment burial for optimal growth
'V_ver' : 0., # [m/year] Vertical growth potential
'germinate' : 0., # [1/year] Possibility of germination per year
'lateral' : 0., # [1/year] Posibility of lateral expension per year
'veg_gamma' : 1., # [-] Constant on influence of sediment burial
'veg_sigma' : 0., # [-] Sigma in gaussian distrubtion of vegetation cover filter
'sedimentinput' : 0., # [-] Constant boundary sediment influx (only used in solve_pieter)
'scheme' : 'euler_backward', # Name of numerical scheme (euler_forward, euler_backward or crank_nicolson)
'solver' : 'trunk', # Name of the solver (trunk, pieter, steadystate,steadystatepieter)
'boundary_lateral' : 'circular', # Name of lateral boundary conditions (circular, constant ==noflux)
'boundary_offshore' : 'constant', # Name of offshore boundary conditions (flux, constant, uniform, gradient)
'boundary_onshore' : 'gradient', # Name of onshore boundary conditions (flux, constant, uniform, gradient)
'boundary_gw' : 'no_flow', # Landward groundwater boundary, dGw/dx = 0 (or 'static')
'method_moist_threshold' : 'belly_johnson', # Name of method to compute wind velocity threshold based on soil moisture content
'method_moist_process' : 'infiltration', # Name of method to compute soil moisture content(infiltration or surface_moisture)
'offshore_flux' : 0., # [-] Factor to determine offshore boundary flux as a function of Q0 (= 1 for saturated flux , = 0 for noflux)
'constant_offshore_flux' : 0., # [kg/m/s] Constant input flux at offshore boundary
'onshore_flux' : 0., # [-] Factor to determine onshore boundary flux as a function of Q0 (= 1 for saturated flux , = 0 for noflux)
'constant_onshore_flux' : 0., # [kg/m/s] Constant input flux at offshore boundary
'lateral_flux' : 0., # [-] Factor to determine lateral boundary flux as a function of Q0 (= 1 for saturated flux , = 0 for noflux)
'method_transport' : 'bagnold', # Name of method to compute equilibrium sediment transport rate
'method_roughness' : 'constant', # Name of method to compute the roughness height z0, note that here the z0 = k, which does not follow the definition of Nikuradse where z0 = k/30.
'method_grainspeed' : 'windspeed', # Name of method to assume/compute grainspeed (windspeed, duran, constant)
'method_shear' : 'fft', # Name of method to compute topographic effects on wind shear stress (fft, quasi2d, duna2d (experimental))
'max_error' : 1e-6, # [-] Maximum error at which to quit iterative solution in implicit numerical schemes
'max_iter' : 1000, # [-] Maximum number of iterations at which to quit iterative solution in implicit numerical schemes
'max_iter_ava' : 1000, # [-] Maximum number of iterations at which to quit iterative solution in avalanching calculation
'refdate' : '2020-01-01 00:00', # [-] Reference datetime in netCDF output
'callback' : None, # Reference to callback function (e.g. example/callback.py':callback)
'wind_convention' : 'nautical', # Convention used for the wind direction in the input files (cartesian or nautical)
'alfa' : 0, # [deg] Real-world grid cell orientation wrt the North (clockwise)
'dune_toe_elevation' : 3, # Choose dune toe elevation, only used in the PH12 dune erosion solver
'beach_slope' : 0.1, # Define the beach slope, only used in the PH12 dune erosion solver
'veg_min_elevation' : -10., # Minimum elevation (m) where vegetation can grow; default -10 disables restriction (allows vegetation everywhere). Set to a higher value to enforce a minimum elevation for vegetation growth.
'vegshear_type' : 'raupach', # Choose the Raupach grid based solver (1D or 2D) or the Okin approach (1D only)
'okin_c1_veg' : 0.48, #x/h spatial reduction factor in Okin model for use with vegetation
'okin_c1_fence' : 0.48, #x/h spatial reduction factor in Okin model for use with sand fence module
'okin_initialred_veg' : 0.32, #initial shear reduction factor in Okin model for use with vegetation
'okin_initialred_fence' : 0.32, #initial shear reduction factor in Okin model for use with sand fence module
'veggrowth_type' : 'orig', #'orig', 'duranmoore14'
'rhoveg_max' : 0.5, #maximum vegetation density, only used in duran and moore 14 formulation
't_veg' : 3, #time scale of vegetation growth (days), only used in duran and moore 14 formulation
'v_gam' : 1, # only used in duran and moore 14 formulation
}
REQUIRED_CONFIG = ['nx', 'ny']
Processes, Thresholds, Methods, and Boundary Conditions
Processes in AeoLiS can be activated (T) or deactivated (F) via the configuration file. Several thresholds (th_...) can be set seperately and methods keywords (method_...) determine which methods are used to compute these processes. Seperate keywords set which boundary conditions are applied.
More detailed descriptions of these processes can be found in the model description. To see exactly where specific process, threshold, and method flags are utilized in the code, an easy approach is to use the search bar on the main page of the AeoLiS GitHub repository.
List overview
process_wind: Interpolates wind speed and direction to each time step. The model cannot run without this flag.process_threshold: Computes the threshold shear velocity based on various processes. (Note: This is overridden if a static threshold file is provided).process_transport: Computes the sediment transport rates.process_bedupdate: Allows the bed level to change based on calculated erosion and deposition.process_shear: Computes the shear stress perturbation caused by topography (needed for landform simulations).process_separation: Computes flow separation bubbles over steep slopes (requiresprocess_shear).process_avalanche: Computes avalanching when bed slopes exceed a critical static angle.process_tide: Computes water level elevations, determining whether cells becomes submerged.process_wave: Computes wave heights across the domain, including $Hs_{mix}$, required for the mixing of sediment.process_runup: Computes the runup extent using the Stockdon equation (requiresprocess_waveandprocess_tide).process_moist: Computes soil moisture content via different surface moisture methods. (Related flags includeprocess_groundwater,process_seepage_face, andprocess_scanning).process_mixtoplayer: Mixes sediment fractions over several layers down to the Depth of Disturbance (requiresprocess_wave).process_wet_bed_reset: Resets the bed to the original bathymetry if submerged. Execution depends on Total Water Level (TWL), requiringprocess_runup,process_wave, andprocess_tide.process_vegetation: Computes vegetation growth and shear stress reduction.process_fences: Computes shear velocity reduction based on user-provided fence characteristics (Okin model).process_dune_erosion: Computes dune erosion and triggers avalanching when water levels impact the dunes.process_meteo: (Placeholder flag; currently has no functionality).
Thresholds
th_grainsize: Computes threshold velocity based on grain size (base value,uth0).th_moisture: Modifies threshold velocity based on moisture content (requiresprocess_moist).th_sheltering: Modifies threshold based on sheltering by roughness elements (coarser sediment fractions).th_nelayer: Modifies threshold based on the presence of a non-erodible layer.(Note: ``th_bedslope``, ``th_drylayer``, ``th_salt``, and ``th_humidity`` are currently not fully implemented or rarely used).
Methods
Note: More guidance on the decision for ``method_grainspeed``, ``method_shear``, and ``solver`` is given in the “Guidance on advection, shear and grainspeed solvers” section later on.
method_transport: Defines the transport equation (Options:bagnold,bagnold_gs,kawamura,lettau,dk,sauermann,vanrijn_strypsteen).method_grainspeed: Defines the speed of sediment in the air (Options:windspeed,duran,duran_uniform,duran_full).method_shear: Computes topographic effects on wind shear stress (Options:fft,1Dstacks).method_roughness: Computes the roughness height (Options:constant,constant_nikuradse,vanrijn_strypsteen).method_vegetation: Defines the vegetation formulation (Options:duran,grass). The newest vegetation implementation is called throughgrass.method_moist_process: Computes soil moisture content (Options:infiltration,surface_moisture).method_moist_threshold: Computes wind velocity threshold based on soil moisture (Options:belly_johnson).solver: Defines the numerical advection solver (Options:steadystate,euler_backward,euler_forward).
Boundary Conditions
boundary_offshore,boundary_onshore,boundary_lateral: Defines the method for handling sediment concentrations at the domain edges. Options include:constant: Applies a zero-gradient boundary condition (incoming concentration equals the adjacent internal cell).flux: Applies a user-defined incoming sediment flux based on the equilibrium concentration.circular: Applies periodic boundaries where sediment leaving one side re-enters the opposite side (Note: offshore and onshore boundaries must both be set tocirculartogether).
offshore_flux,onshore_flux,lateral_flux: Defines the incoming sediment concentration as a fraction of the equilibrium concentration (Cu) when the respective boundary is set toflux(e.g.,1.0for fully saturated incoming wind,0.0for clean air with no sediment).
Typical Configurations
When building an AeoLiS model, it helps to know which processes and methods are typically required for different environments. The table below provides a general baseline for five common simulation types.
Parameter / Flag |
Wind-only [1] |
Barchan [2] |
Parabolic [3] |
Beach [4] |
Blowout [5] |
|---|---|---|---|---|---|
Processes |
|||||
|
✅ |
✅ |
✅ |
✅ |
✅ |
|
✅ |
✅ |
✅ |
✅ |
✅ |
|
✅ |
✅ |
✅ |
✅ |
✅ |
|
✅ |
✅ |
✅ |
✅ |
✅ |
|
❌ |
✅ |
✅ |
✅ |
✅ |
|
❌ |
✅ |
✅ |
❌ |
❌ |
|
❌ |
✅ |
✅ |
✅ |
✅ |
|
❌ |
❌ |
✅ |
✅ |
✅ |
|
❌ |
❌ |
❌ |
✅ |
✅ |
|
❌ |
❌ |
❌ |
✅ |
✅ |
|
❌ |
❌ |
❌ |
✅ |
✅ |
|
❌ |
❌ |
❌ |
✅ |
✅ |
|
❌ |
❌ |
❌ |
✅ |
✅ |
|
❌ |
❌ |
❌ |
✅ |
✅ |
|
❌ |
❌ |
❌ |
✅ (maybe) |
❌ |
Thresholds |
|||||
|
✅ |
✅ |
✅ |
✅ |
✅ |
|
❌ |
❌ |
❌ |
✅ |
✅ |
|
❌ |
❌ |
❌ |
✅ |
✅ |
|
❌ |
✅ |
✅ |
❌ |
❌ |
Methods |
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
Boundary Conditions |
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Model output
AeoLiS writes its simulation results to a NetCDF4 file (by default named aeolis.nc). The outputed variables (output_vars) and frequency (output_times) of this file are configured via the configuration file. You can request any of the variables listed in the Model state/output section, while x and y are automatically included.
While most output variables have 2D spatial dimensions combined with a time dimension (time, ny, nx), some contain additional dimensions. For example, sediment mass (mass) is calculated per grid cell, bed layer, and sediment fraction. Outputting such a five-dimensional variable can result in large file sizes. To reduce this, AeoLiS offers masstop as output, which provides the sediment distribution only for the active top layer.
By default, the output represents the model state at the exact moment defined by the output interval. If the internal time step (dt) is smaller than this interval, intermediate calculations are not saved. To evaluate variable behavior between output intervals, statistical summaries can be requested by appending suffixes directly to the variable names in output_vars. Available suffixes include _avg (average), _sum (cumulative sum), _var (variance), _min (minimum), and _max (maximum) over the output interval.
Note
Known Issue: Because the model parses the underscore (_) to identify these statistical requests, variables that contain underscores in their base names will cause parsing conflicts.
Example Python script
Below is a Python script showing how to read the NetCDF output data en create a single plot.
import netCDF4 as nc
import matplotlib.pyplot as plt
# 1. Open the NetCDF output file
ncfile = 'aeolis.nc'
ds = nc.Dataset(ncfile, 'r')
# 2. Load spatial coordinates and bed level
x = ds.variables['x'][:, :]
y = ds.variables['y'][:, :]
zb_final = ds.variables['zb'][-1, :, :] # (time, y, x)
ds.close()
# 3. Plot the final bed level
fig, ax = plt.subplots(figsize=(8, 6))
pc = ax.pcolormesh(x, y, zb_final, cmap='viridis', shading='auto')
ax.set_aspect('equal')
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_title('Final Bed Level (zb)')
fig.colorbar(pc, ax=ax, label='Bed level (m)')
plt.tight_layout()
plt.show()
Example Python script for animation
Below is a second example demonstrating how to animate time-series output of the cross-shore (ustars) and alongshore (ustarn) shear velocity components.
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# 1. Open the NetCDF output file and load variables
ncfile = 'aeolis.nc'
ds = nc.Dataset(ncfile, 'r')
time = ds.variables['time'][:]
x = ds.variables['x'][:, :]
y = ds.variables['y'][:, :]
ustar = ds.variables['ustar'][:] # (time, y, x)
ustars = ds.variables['ustars'][:] # (time, y, x)
ustarn = ds.variables['ustarn'][:] # (time, y, x)
ds.close()
# 2. Set up the figure layout
fig, ax = plt.subplots(figsize=(8, 6))
title = ax.set_title(f'Shear Velocity at t = {time[0]:.0f} s')
ax.set_aspect('equal')
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
# Initial background mesh, colorbar and overlaying quiver plot (vectors)
pc = ax.pcolormesh(x, y, ustar[0, :, :], cmap='YlOrRd', shading='auto', vmin=0, vmax=np.max(ustar))
fig.colorbar(pc, ax=ax, label='Shear velocity magnitude (m/s)')
Q = ax.quiver(x, y, ustars[0, :, :], ustarn[0, :, :], color='black')
# 3. Define the update function for the animation
def update(frame):
pc.set_array(ustar[frame, :, :].ravel()) # update basemap
Q.set_UVC(ustars[frame, :, :], ustarn[frame, :, :]) # update quiver
title.set_text(f'Shear Velocity at t = {time[frame]:.0f} s') # update title
return pc, Q, title
# 4. Create the animation
ani = animation.FuncAnimation(fig, update, frames=len(time), blit=False)
ani.save('ustar_animation.mp4', writer='ffmpeg', fps=10)
Model state/output
The AeoLiS model state is described by a collection of spatial grid variables with at least one value per horizontal grid cell. Specific model state variables can also be subdivided over bed composition layers and/or grain size fractions. All model state variables can be part of the model netCDF4 output. The current model state variables are listed below.
INITIAL_STATE = {
('ny', 'nx') : (
'uw', # [m/s] Wind velocity
'uws', # [m/s] Component of wind velocity in x-direction
'uwn', # [m/s] Component of wind velocity in y-direction
'tau', # [N/m^2] Wind shear stress
'taus', # [N/m^2] Component of wind shear stress in x-direction
'taun', # [N/m^2] Component of wind shear stress in y-direction
'tau0', # [N/m^2] Wind shear stress over a flat bed
'taus0', # [N/m^2] Component of wind shear stress in x-direction over a flat bed
'taun0', # [N/m^2] Component of wind shear stress in y-direction over a flat bed
'taus_u', # [N/m^2] Saved direction of wind shear stress in x-direction
'taun_u', # [N/m^2] Saved direction of wind shear stress in y-direction
'dtaus', # [-] Component of the wind shear perturbation in x-direction
'dtaun', # [-] Component of the wind shear perturbation in y-direction
'ustar', # [m/s] Wind shear velocity
'ustars', # [m/s] Component of wind shear velocity in x-direction
'ustarn', # [m/s] Component of wind shear velocity in y-direction
'ustar0', # [m/s] Wind shear velocity over a flat bed
'ustars0', # [m/s] Component of wind shear velocity in x-direction over a flat bed
'ustarn0', # [m/s] Component of wind shear velocity in y-direction over a flat bed
'udir', # [rad] Wind direction
'zs', # [m] Water level above reference (or equal to zb if zb > zs)
'SWL', # [m] Still water level above reference
'Hs', # [m] Wave height
'Hsmix', # [m] Wave height for mixing (including setup, TWL)
'Tp', # [s] Wave period for wave runup calculations
'zne', # [m] Non-erodible layer
),
}
MODEL_STATE = {
('ny', 'nx') : (
'x', # [m] Real-world x-coordinate of grid cell center
'y', # [m] Real-world y-coordinate of grid cell center
'ds', # [m] Real-world grid cell size in x-direction
'dn', # [m] Real-world grid cell size in y-direction
'dsdn', # [m^2] Real-world grid cell surface area
'dsdni', # [m^-2] Inverse of real-world grid cell surface area
# 'alfa', # [rad] Real-world grid cell orientation #Sierd_comm in later releases this needs a revision
'zb', # [m] Bed level above reference
'zs', # [m] Water level above reference
'zne', # [m] Height above reference of the non-erodible layer
'zb0', # [m] Initial bed level above reference
'zdry', # [m]
'dzdry', # [m]
'dzb', # [m/dt] Bed level change per time step (computed after avalanching!)
'dzbyear', # [m/yr] Bed level change translated to m/y
'dzbavg', # [m/year] Bed level change averaged over collected time steps
'S', # [-] Level of saturation
'moist', # [-] Moisture content (volumetric)
'moist_swr', # [-] Moisture content soil water retention relationship (volumetric)
'h_delta', # [-] Suction at reversal between wetting/drying conditions
'gw', # [m] Groundwater level above reference
'gw_prev', # [m] Groundwater level above reference in previous timestep
'wetting', # [bool] Flag indicating wetting or drying of soil profile
'scan_w', # [bool] Flag indicating that the moisture is calculated on the wetting scanning curve
'scan_d', # [bool] Flag indicating that the moisture is calculated on the drying scanning curve
'scan_w_moist', # [-] Moisture content (volumetric) computed on the wetting scanning curve
'scan_d_moist', # [-] Moisture content (volumetric) computed on the drying scanning curve
'w_h', # [-] Moisture content (volumetric) computed on the main wetting curve
'd_h', # [-] Moisture content (volumetric) computed on the main drying curve
'w_hdelta', # [-] Moisture content (volumetric) computed on the main wetting curve for hdelta
'd_hdelta', # [-] Moisture content (volumetric) computed on the main drying curve for hdelta
'ustar', # [m/s] Shear velocity by wind
'ustars', # [m/s] Component of shear velocity in x-direction by wind
'ustarn', # [m/s] Component of shear velocity in y-direction by wind
'ustar0', # [m/s] Initial shear velocity (without perturbation)
'zsep', # [m] Z level of polynomial that defines the separation bubble
'hsep', # [m] Height of separation bubble = difference between z-level of zsep and of the bed level zb
'theta_dyn', # [degrees] spatially varying dynamic angle of repose for avalanching
'rhoveg', # [-] Vegetation cover
'drhoveg', # Change in vegetation cover
'hveg', # [m] height of vegetation
'dhveg', # [m] Difference in vegetation height per time step
'dzbveg', # [m] Bed level change used for calculation of vegetation growth
'germinate', # [bool] Newly vegetated due to germination (or establishment)
'lateral', # [bool] Newly vegetated due to lateral propagation
'vegetated', # [bool] Vegetated, determines if vegetation growth or burial is allowed
'vegfac', # Vegetation factor to modify shear stress by according to Raupach 1993
'fence_height', # Fence height
'R', # [m] wave runup
'eta', # [m] wave setup
'sigma_s', # [m] swash
'TWL', # [m] Total Water Level above reference (SWL + Run-up)
'SWL', # [m] Still Water Level above reference
'DSWL', # [m] Dynamic Still water level above reference (SWL + Set-up)
'Rti', # [-] Factor taking into account sheltering by roughness elements
),
('ny','nx','nfractions') : (
'Cu', # [kg/m^2] Equilibrium sediment concentration integrated over saltation height
'Cuf', # [kg/m^2] Equilibrium sediment concentration integrated over saltation height, assuming the fluid shear velocity threshold
'Cu0', # [kg/m^2] Flat bad equilibrium sediment concentration integrated over saltation height
'Ct', # [kg/m^2] Instantaneous sediment concentration integrated over saltation height
'q', # [kg/m/s] Instantaneous sediment flux
'qs', # [kg/m/s] Instantaneous sediment flux in x-direction
'qn', # [kg/m/s] Instantaneous sediment flux in y-direction
'pickup', # [kg/m^2] Sediment entrainment
'w', # [-] Weights of sediment fractions
'w_init', # [-] Initial guess for ``w''
'w_air', # [-] Weights of sediment fractions based on grain size distribution in the air
'w_bed', # [-] Weights of sediment fractions based on grain size distribution in the bed
'uth', # [m/s] Shear velocity threshold
'uthf', # [m/s] Fluid shear velocity threshold
'uth0', # [m/s] Shear velocity threshold based on grainsize only (aerodynamic entrainment)
'u', # [m/s] Mean horizontal saltation velocity in saturated state
'us', # [m/s] Component of the saltation velocity in x-direction
'un', # [m/s] Component of the saltation velocity in y-direction
'usST', # [NEW] [m/s] Component of the saltation velocity in x-direction for SedTRAILS
'unST', # [NEW] [m/s] Component of the saltation velocity in y-direction for SedTRAILS
'u0',
'masstop', # [kg/m^2] Sediment mass in bed toplayer, only stored for output
),
('ny','nx','nlayers') : (
'thlyr', # [m] Bed composition layer thickness
'salt', # [-] Salt content
),
('ny','nx','nlayers','nfractions') : (
'mass', # [kg/m^2] Sediment mass in bed
),
}
Guide on model schematization
Advection solver
The advection equation is the core of sediment transport computations in AeoLiS. It balances the spatial change in sediment transport with the pickup or deposition rate:
Where \(C\) is the actual sediment concentration, \(U_s\) is the sediment velocity, \(x\) is the distance in the transport direction, \(C_{\mathrm{sat}}\) is the equilibrium sediment concentration, and \(T\) is the adaptation time scale.
You can solve this using three different methods (via the solver keyword):
steadystate (default): Highly recommended for most cases. Assumes concentration does not change over time within a single timestep (\(dc/dt = 0\)) and solves rapidly using a finite difference sweeping algorithm.
euler_backward: An implicit, time-varying solver (formerly
pieter_solver). Recommended primarily for legacy purposes or when the steady-state assumption is problematic (e.g., \(dt < 10\) s).euler_forward: A simple explicit solver. Slow due to strict CFL stability conditions. Included mainly for testing or educational purposes.
Shear and grainspeed schematization
There are several options for schematizing spatial shear and sediment velocity. Higher complexity yields more realistic physics but increases computational cost. To illustrate the impact of these choices, we use a consistent demonstration case: A cone-shaped bedform sits on a non-erodible layer, upwind of a vegetated patch designed to capture all incoming sediment. We expect the landform to migrate downwind and evolve into a crescentic barchan dune. We also track the amount of deposition in the vegetated area to monitor mass balance.
Case 0: Windspeed - original method by Hoonhout (2016)
Configuration:
process_shear = F,method_grainspeed = windspeed
Topographic feedback is disabled (\(\nabla \tau = 0\)), meaning the wind blows over the cone as if it were flat. The sediment velocity simply equals the wind speed (\(u_s = u_w\)). As a result, no landform development occurs. Sediment is stripped from the cone and deposited directly into the vegetation. This is the fastest method (20 mins), but should only be used for bulk transport calculations where morphodynamics are irrelevant.
Case 1: Uniform grainspeed
Configuration:
process_shear = F,method_grainspeed = duran_uniform
Sediment in saltation moves significantly slower than the wind. This case introduces the saltation model by Durán (2007) so that \(u_s < u_w\), but keeps the flat-bed assumption for shear velocity. Deposition patterns become more localized, but the landform still does not migrate because there is no topographic steering. It is slightly slower (33 mins) but physically more realistic than Case 0 for static topographies.
Case 2: Topographic steering
Configuration:
process_shear = T,method_grainspeed = duran_uniform
The shear velocity vector is now perturbed by the topography, creating spatial gradients in transport capacity. The sediment velocity \(u_s\) remains uniform. The variation in shear velocity causes the entire landform to migrate towards the vegetation. However, without flow separation or spatial variation in grain speed, it fails to evolve into a crescentic shape. Computational cost increases significantly (2:01 hrs) due to the secondary rotational grid required to calculate shear and increased spatial variability in transport which causes the advection solver to need more iterations.
Case 3: Separation bubble (special)
Configuration:
process_separation = T(Requiresprocess_shear = T)
Steep lee-side slopes cause airflow to separate, creating a zone of recirculation and low shear. Activating the separation bubble prevents sediment transport on the lee side, preserving the steep slip face and allowing the crescentic barchan shape to form (2:10 hrs). Note: In highly complex topographies (like dense vegetation), the bubble may produce undesirable morphodynamics, so use it judiciously.
Case 4: Spatially varying grainspeed
Configuration:
process_shear = T,method_grainspeed = duran
For a complete description of transport around landforms, the sediment velocity \(u_s\) must also vary spatially. This method uses Durán’s analytical approximation, which incorporates slope terms but assumes slopes are relatively gentle to avoid heavy numerical solving. Produces the expected crescentic morphodynamics efficiently (2:03 hrs). This is the recommended configuration for most simulations involving bedform evolution where topographic steering is important.
Case 5: Steeper Slopes (numerical solution)
Configuration:
process_shear = T,method_grainspeed = duran_full
The analytical approximation in Case 4 can overestimate upslope transport on very steep slopes (e.g., \(> 33^\circ\)). This method solves the full sediment velocity equation numerically to account for strong gravitational effects. While differences are subtle in standard dune simulations, this method is crucial for extreme topography like steep blowout cliffs (2:09 hrs). Caution: Dynamic avalanching is not fully coupled yet, so schematize steep slopes carefully.
Summary of Parameter Settings
Case |
Description |
|
|
|
Application |
Time |
|---|---|---|---|---|---|---|
0 |
Flat / Windspeed |
|
|
|
Bulk transport |
0:20 |
1 |
Flat / Grainspeed |
|
|
|
Static topography |
0:33 |
2 |
Topo Steering |
|
|
|
Veg-dominated dunes |
2:01 |
3 |
Separation Bubble |
|
|
|
Moderate morphology |
2:10 |
4 |
Varying $u_s$ (Approx) |
|
|
|
Landform evolution |
2:03 |
5 |
Varying $u_s$ (Full) |
|
|
|
Extreme topography |
2:09 |