Simple wave sediment modelling

Note

Using a wave propagation model based on linear wave theory, we investigate the cumulative effect of a given wave forcing condition on erosional and depositional trend across a regional coastal system.

Loading required modules for this exercise. The main module here is wavesed a piece of code that was developed for this course and that can be found on GitHub

!cd wavesed/; f2py -c -m ocean ocean.f90    
running build
running config_cc
unifing config_cc, config, build_clib, build_ext, build commands --compiler options
running config_fc
unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options
running build_src
build_src
building extension "ocean" sources
f2py options: []
f2py:> /tmp/tmp09p_7588/src.linux-x86_64-3.6/oceanmodule.c
creating /tmp/tmp09p_7588/src.linux-x86_64-3.6
Reading fortran codes...
	Reading file 'ocean.f90' (format:free)
Post-processing...
	Block: ocean
			Block: ocean
				Block: airymodel
				Block: transport
				Block: diffusion
Post-processing (stage 2)...
	Block: ocean
		Block: unknown_interface
			Block: ocean
				Block: airymodel
				Block: transport
				Block: diffusion
Building modules...
	Building module "ocean"...
		Constructing F90 module support for "ocean"...
		  Variables: pi onpi pi2 pion2 onpi2 grav
			Constructing wrapper function "ocean.airymodel"...
			  c,l,travel,waveh = airymodel(dx,dd,h0,depth,src,inland,shadow,[numrow,numcol])
			Constructing wrapper function "ocean.transport"...
			  dz,dist = transport(iter,depth,hent,transx,transy,[numrow,numcol])
			Constructing wrapper function "ocean.diffusion"...
			  depo = diffusion(oelev,dz,coeff,maxth,tstep,nstep,[numrow,numcol])
	Wrote C/API module "ocean" to file "/tmp/tmp09p_7588/src.linux-x86_64-3.6/oceanmodule.c"
	Fortran 90 wrappers are saved to "/tmp/tmp09p_7588/src.linux-x86_64-3.6/ocean-f2pywrappers2.f90"
  adding '/tmp/tmp09p_7588/src.linux-x86_64-3.6/fortranobject.c' to sources.
  adding '/tmp/tmp09p_7588/src.linux-x86_64-3.6' to include_dirs.
copying /usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmp09p_7588/src.linux-x86_64-3.6
copying /usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmp09p_7588/src.linux-x86_64-3.6
  adding '/tmp/tmp09p_7588/src.linux-x86_64-3.6/ocean-f2pywrappers2.f90' to sources.
build_src: building npy-pkg config files
running build_ext
customize UnixCCompiler
C compiler: /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-cc -DNDEBUG -fwrapv -O2 -Wall -Wstrict-prototypes -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/coast/include -DNDEBUG -D_FORTIFY_SOURCE=2 -O2 -isystem /usr/share/miniconda/envs/coast/include -fPIC

creating /tmp/tmpcnr870fu/tmp
creating /tmp/tmpcnr870fu/tmp/tmpcnr870fu
compile options: '-MMD -MF /tmp/tmpcnr870fu/file.c.d -c'
x86_64-conda-linux-gnu-cc: /tmp/tmpcnr870fu/file.c
customize UnixCCompiler using build_ext
get_default_fcompiler: matching types: '['gnu95', 'intel', 'lahey', 'pg', 'nv', 'absoft', 'nag', 'vast', 'compaq', 'intele', 'intelem', 'gnu', 'g95', 'pathf95', 'nagfor']'
customize Gnu95FCompiler
Found executable /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-gfortran
Found executable /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-ld
Found executable /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-ar
Found executable /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-ranlib
customize Gnu95FCompiler
customize Gnu95FCompiler using build_ext
building 'ocean' extension
compiling C sources
C compiler: /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-cc -DNDEBUG -fwrapv -O2 -Wall -Wstrict-prototypes -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/coast/include -DNDEBUG -D_FORTIFY_SOURCE=2 -O2 -isystem /usr/share/miniconda/envs/coast/include -fPIC

creating /tmp/tmp09p_7588/tmp
creating /tmp/tmp09p_7588/tmp/tmp09p_7588
creating /tmp/tmp09p_7588/tmp/tmp09p_7588/src.linux-x86_64-3.6
compile options: '-I/tmp/tmp09p_7588/src.linux-x86_64-3.6 -I/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/include -I/usr/share/miniconda/envs/coast/include/python3.6m -c'
x86_64-conda-linux-gnu-cc: /tmp/tmp09p_7588/src.linux-x86_64-3.6/oceanmodule.c
x86_64-conda-linux-gnu-cc: /tmp/tmp09p_7588/src.linux-x86_64-3.6/fortranobject.c
In file included from /usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1822,
                 from /usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:12,
                 from /usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /tmp/tmp09p_7588/src.linux-x86_64-3.6/fortranobject.h:13,
                 from /tmp/tmp09p_7588/src.linux-x86_64-3.6/fortranobject.c:2:
/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:17:2: warning: #warning "Using deprecated NumPy API, disable it with " "#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
   17 | #warning "Using deprecated NumPy API, disable it with " \
      |  ^~~~~~~
In file included from /usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1822,
                 from /usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:12,
                 from /usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /tmp/tmp09p_7588/src.linux-x86_64-3.6/fortranobject.h:13,
                 from /tmp/tmp09p_7588/src.linux-x86_64-3.6/oceanmodule.c:16:
/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:17:2: warning: #warning "Using deprecated NumPy API, disable it with " "#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
   17 | #warning "Using deprecated NumPy API, disable it with " \
      |  ^~~~~~~
/tmp/tmp09p_7588/src.linux-x86_64-3.6/oceanmodule.c:109:12: warning: 'f2py_size' defined but not used [-Wunused-function]
  109 | static int f2py_size(PyArrayObject* var, ...)
      |            ^~~~~~~~~
compiling Fortran 90 module sources
Fortran f77 compiler: /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -fopenmp -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/coast/include -O3 -funroll-loops
Fortran f90 compiler: /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-gfortran -Wall -g -fno-second-underscore -fPIC -fopenmp -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/coast/include -O3 -funroll-loops
Fortran fix compiler: /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-gfortran -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -fPIC -fopenmp -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/coast/include -O3 -funroll-loops
compile options: '-I/tmp/tmp09p_7588/src.linux-x86_64-3.6 -I/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/include -I/usr/share/miniconda/envs/coast/include/python3.6m -c'
extra options: '-J/tmp/tmp09p_7588/ -I/tmp/tmp09p_7588/'
x86_64-conda-linux-gnu-gfortran:f90: ocean.f90
compiling Fortran sources
Fortran f77 compiler: /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -fopenmp -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/coast/include -O3 -funroll-loops
Fortran f90 compiler: /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-gfortran -Wall -g -fno-second-underscore -fPIC -fopenmp -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/coast/include -O3 -funroll-loops
Fortran fix compiler: /usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-gfortran -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -fPIC -fopenmp -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/coast/include -O3 -funroll-loops
compile options: '-I/tmp/tmp09p_7588/src.linux-x86_64-3.6 -I/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/include -I/usr/share/miniconda/envs/coast/include/python3.6m -c'
extra options: '-J/tmp/tmp09p_7588/ -I/tmp/tmp09p_7588/'
x86_64-conda-linux-gnu-gfortran:f90: /tmp/tmp09p_7588/src.linux-x86_64-3.6/ocean-f2pywrappers2.f90
/usr/share/miniconda/envs/coast/bin/x86_64-conda-linux-gnu-gfortran -Wall -g -Wall -g -shared -Wl,-O2 -Wl,--sort-common -Wl,--as-needed -Wl,-z,relro -Wl,-z,now -Wl,--disable-new-dtags -Wl,--gc-sections -Wl,-rpath,/usr/share/miniconda/envs/coast/lib -Wl,-rpath-link,/usr/share/miniconda/envs/coast/lib -L/usr/share/miniconda/envs/coast/lib /tmp/tmp09p_7588/tmp/tmp09p_7588/src.linux-x86_64-3.6/oceanmodule.o /tmp/tmp09p_7588/tmp/tmp09p_7588/src.linux-x86_64-3.6/fortranobject.o /tmp/tmp09p_7588/ocean.o /tmp/tmp09p_7588/tmp/tmp09p_7588/src.linux-x86_64-3.6/ocean-f2pywrappers2.o -L/usr/share/miniconda/envs/coast/bin/../x86_64-conda-linux-gnu/sysroot/lib/../lib -L/usr/share/miniconda/envs/coast/bin/../x86_64-conda-linux-gnu/sysroot/lib/../lib -lgfortran -o ./ocean.cpython-36m-x86_64-linux-gnu.so
Removing build directory /tmp/tmp09p_7588
import time
import numpy as np
import wavesed.wave as ocean

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

Model set-up

The model requires several parameters to be defined before being run. The most important one being the initial bathymetry. For this example, we will use the region around Fraser Island at a resolution of 1000 m. The digital elevation model (DEM) comes from the Project 3DGBR.

fraser

Fig. 25 Fraser Island from Google Earth.

The file is called FraserDEM.csv and is located in the data folder. It is made of 3 columns containing the X, Y, and Z coordinates.

Model domain / grid parameters

It is possible to change the initial resolution by using the rfac parameter, as an example:

  • rfac = 1 corresponds to a model at 1000 m resolution (similar to the initial DEM)

  • rfac = 2 transforms the dataset to 2000 m resolution

# Bathymetric filename
bfile = '../pracenv/dataset/FraserDEM.csv'
# Resolution factor
rfac = 1 

Definition of wave parameters

Now we need to define the forcing parameters for our simulation:

  • H0 the deep water wave height in [m],

  • dir the deep water wave source direction at the boundary in degrees (angle is defined counterclock wise from horizontal axis)

  • wbase the maximum depth for wave influence on seabed [m]

  • slvl the sea-level position in [m]

# Wave height
H0 = 2
# Wave direction at boundary 
dir = 0
# Maximum depth for wave influence (m)
wbase = 10
# Sea level position (m)
slvl = 0.
# Output figure name
outname = 'wdir0'

Definition of sediment parameters

To estimate wave-induced sediment transport, this simple model only accounts for a single grain size that needs to be defined (d50).

There are also additional parameters that can be set to perform longshore and diffusive transport but we will use the default values for this example…

# Mean grain size diameter in m
d50 = 0.0001

Model initialisation function

Now that all the parameters are defined in our python environment, we can initialise the simulation by calling the wavesed module:

#help(ocean.wave.__init__)
wavesed = ocean.wave(filename=bfile, wavebase=wbase, resfac=rfac, dia=d50)

As you can see from the command above there are still some parameters that have not been passed to the module:

  • slvl

  • dir

  • H0

These parameters are sent separately so that multiple climatic conditions could be simulate more easily.

First we define the sea-level position slvl.

t0 = time.clock()
wavesed.findland(slvl)
print('Definition of shoreline position took (s):',time.clock()-t0)
Definition of shoreline position took (s): 0.0030769999999999964

Initialising wave boundary direction dir is then done with the wavesource function:

wdir = wavesed.wavesource(dir)

Waves evolution

The waves are then transformed using the linear wave theory from deep (set-up conditions) to shallow water assuming shore-parallel depth contours.

The orientation of wave fronts is determined by wave celerity and refraction due to depth variations.

Travel time in the domain is calculated from Huygen’s principle (using an order \(\sqrt{5}\) approximation).

Assuming no refraction or loss of energy due to bottom friction, wave power \(P\) is conserved from deep to shallow water.

#help(wavesed.cmptwaves)
t0 = time.clock()
wavesed.cmptwaves(src=wdir, h0=H0, sigma=1.)
print('Wave evolution computation took (s): ',time.clock()-t0)
Wave evolution computation took (s):  0.6260160000000001

Plotting simulated wave characteristics

A series of graph can be visualised to show the following characteristics:

  • bathy: coarse bathymetric map

  • wlength: coarse wave lenght map

  • travel: coarse wave travel time map

  • wcelerity: coarse wave celerity map

  • wheight: coarse wave height map

  • wpower: coarse wave power map

  • ubot: coarse wave induced bottom velocity map

  • shear: coarse wave induced shear stress map

#help(wavesed.plotData)
size = (6,6)
# i1 = 0 
# i2 = -1
# j1 = 0
# j2 = -1

# Zooming to a specific region
i1 = 80 
i2 = 250 
j1 = 0
j2 = 170 

wavesed.plotData(data='bathy', figsize=size, vmin=-100, vmax=100, 
                 fontsize=10, imin=i1, imax=i2, jmin=j1, jmax=j2)

wavesed.plotData(data='travel', figsize=size, tstep=100, vmin=0, vmax=5.e4, 
                 fontsize=10, imin=i1, imax=i2, jmin=j1, jmax=j2)

wavesed.plotData(data='wcelerity', figsize=size, vmin=0, vmax=15, 
                 fontsize=10, stream=3, imin=i1, imax=i2, jmin=j1, jmax=j2)

wavesed.plotData(data='ubot', figsize=size, vmin=0, vmax=2, 
                 fontsize=10, imin=i1, imax=i2, jmin=j1, jmax=j2)

wavesed.plotData(data='shear', figsize=size, vmin=-10, vmax=10, 
                  fontsize=10, imin=i1, imax=i2, jmin=j1, jmax=j2)
../_images/sedtransport_24_0.svg../_images/sedtransport_24_1.svg../_images/sedtransport_24_2.svg../_images/sedtransport_24_3.svg../_images/sedtransport_24_4.svg

Sediment entrainment, transport & deposition

Sediment entrainment relates to wave induced shear stress.

The transport is computed according to both wave direction and longshore transport. Deposition is dependent of shear stress and diffusion.

#help(wavesed.cmptsed)
t0 = time.clock()
wavesed.cmptsed(sigma=1.)
print('Sediment erosion/deposition computation took (s): ',time.clock()-t0)
Sediment erosion/deposition computation took (s):  0.5779190000000014

Plotting morphological changes

Using the same plotting function (plotData), we can now plot the following metrics:

  • ent: coarse sediment entrainment map

  • dz: coarse erosion/deposition map

#help(wavesed.plotData)
size = (6,6)
# i1 = 0 
# i2 = -1
# j1 = 0
# j2 = -1

# Zooming to a specific region
i1 = 80 
i2 = 250 
j1 = 0
j2 = 170

wavesed.plotData(data='bathy', figsize=size, vmin=-100, vmax=100, 
                 fontsize=10, stream=0, imin=i1, imax=i2, jmin=j1, jmax=j2)

wavesed.plotData(data='ent', figsize=size, vmin=0., vmax=5., 
                 fontsize=10, stream=0, imin=i1, imax=i2, jmin=j1, jmax=j2)

wavesed.plotData(data='dz', figsize=size, vmin=-1., vmax=1., 
                 fontsize=10, stream=0, imin=i1, imax=i2, jmin=j1, jmax=j2) #, save=outname)
../_images/sedtransport_30_0.svg../_images/sedtransport_30_1.svg../_images/sedtransport_30_2.svg

Saving wave/sedimentation data

#waveparams.outputCSV(filename='erodep.csv', seddata=1)

Note

  1. Change the initial wave direction conditions to 90 and then to 270. Evaluate the impact on the morphological patterns.

  2. Change the resolution factor to see how it impacts on erosion/deposition resolution.