- Email: [email protected]
Stretched Eulerian coordinate model of coastal sediment transport
ARTICLE IN PRESS Computers & Geosciences 35 (2009) 1167–1176
Contents lists available at ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
Stretched Eulerian coordinate model of coastal sediment transport N. Margvelashvili CSIRO Marine and Atmospheric Research, GPO Box 1538, Hobart, Tasmania 7001, Australia
a r t i c l e in fo
abstract
Article history: Received 16 January 2008 Accepted 28 March 2008
A one-dimensional vertical model of sediment and dissolved material transport in a coupled, vertically resolved benthic-pelagic system has been developed. An apparently unique feature of the model is that the governing equations are first formulated in a time-varying Eulerian frame and then solved numerically, which provides a consistent approach for incorporating resuspension/deposition, swelling/consolidation, and biodiffusion processes within a single modelling framework. Numerical experiments were carried out with varying grids to test diffusive properties of the model. A stretched numerical grid that maintains high vertical resolution of the top sediments and constrains numerical diffusion in deep layers is an integral part of the model. The developed model provides a test bed for testing scientific hypothesis pertaining to coastal sediment transport. The model can also be used to support biogeochemical studies in a coupled benthic-pelagic system. Crown Copyright & 2008 Published by Elsevier Ltd. All rights reserved.
Keywords: Coastal Sediment Numerical Model Benthic Pelagic
1. Introduction Vertical profiles of physical, chemical and biological properties in coastal water and sediments often exhibit strong spatial gradients near the sediment–water interface (Dyer, 1989; Boudreau, 1997). The location of the interface varies with time due to resuspension and deposition or swelling/consolidation of the seabed, and the displacement of this boundary may exceed typical scales of spatial variability of the vertical profiles. The dynamic nature and high spatial variability of the nearbottom fields provide a challenge for adequate description of benthic fluxes across the water and sediment layers, which is essential for better understanding and efficient management of coastal ecosystems. To address this challenge, we developed a modelling framework for coastal material transport that accommodates key sediment processes and reconciles strong vertical gradients with the dynamic nature of coastal sediments.
Tel.: +61 3 62325142; fax: +61 3 62325123.
E-mail address: [email protected]
Integrated coastal models typically include ecological reaction models and physical transport models, operating in a coupled benthic-pelagic system (Skogen et al., 1995; Wainright and Hopkinson, 1997; Luff and Moll, 2004). Transport models in sediments typically employ either a time-invariant Eulerian frame or Lagrangian coordinates. Lagrangian coordinates can provide high vertical resolution of the uppermost sediments during the simulation (Gibson et al., 1967; Govindaraju et al., 1999). However these coordinates are not well suited for simulating mixtures of sand, fines and dissolved materials (which is typically the case with natural sediments), since in a Lagrangian frame the coordinates are fixed to a given material parcel which should always embrace the same solid or dissolved particles, a particular choice being dependent on the simulation objectives (Narasimhan, 1992). In Eulerian coordinates, having all grid layers fixed would imply constant thickness of sediments. To adjust a numerical grid to the time-varying thickness, in Eulerian coordinates either new sediment layers must be added on top of existing layers or old layers must be diluted during the course of simulation. To maintain a high speed in the simulation, the resolution of top sediments must vary
0098-3004/$ - see front matter Crown Copyright & 2008 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2008.03.016
ARTICLE IN PRESS 1168
N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
dramatically with time, compromising the accuracy of the calculations. An alternative approach for simulating material transport in sediments is to solve transport equations using Eulerian grids, where all grid layers vary with time (Harris and Wiberg, 2001; Reed et al., 1999). Since the time-varying grid layers introduce large numerical diffusion in sediments, these models simulate suspended sediment transport over a two-layered seabed assuming a fully mixed environment under the top sediments. The model described in this paper extends the latter formulation to multilayered sediments. It employs timevarying Eulerian grids plus a user-specified dynamics for the grid layers. To maintain high resolution in the top sediments without sacrificing computational efficiency, the model utilises varying time scales of the grid adjustment to the sediment thickness. The developed model conserves masses of solid, liquid and dissolved phases, and provides a consistent approach for incorporating resuspension/deposition, swelling/consolidation, and biodiffusion processes within a single modelling framework.
2. Transport equations Here and throughout the paper, solid and liquid phases are assumed to be incompressible. The vertical axis of the coordinate frame is directed upward. Sediments and dissolved tracers are represented by vertical profiles of the sediment concentration Cs(z, t) (kg/(m3 mixture)), and the concentration of the dissolved tracer Cd(z, t) (kg/(m3 water)). The superscripts (s), (w), and (d) denote the solid phase, liquid phase (water), and dissolved phase, respectively. The subscripts (i) and (j) are introduced whenever sediments or dissolved traces are represented by a number of classes. All classes of particles in sediments have the same vertical velocity associated with the sediment consolidation UisUcs(z, t) and the vertical velocity of all dissolved tracers is equal to the water velocity UidUw(z, t). Sediments and water column are represented by discrete, time-varying layers of thickness Zk ¼ z(k+)z(k), with the subscript k indicating the layer number, and (k+) and (k) indexes indicating the top and bottom interfaces between layers, respectively (Fig. 1).
2.1. Mass conservation A general mass balance equation for sediments (C ¼ Cs) or dissolved tracer (C ¼ Cd), in a stretched Eulerian frame z ¼ z(t), can be expressed as follows (Appendix A): @V g @ @ @ ðC jÞ ¼ C jC . (1) j U˜ n @t @z @z @z Here Vg is the velocity of the time-varying coordinate levels; U˜ ¼ ðU V g Þ is the velocity of either a particulate or dissolved tracer expressed in a time-varying frame z ¼ z(t); U is the tracer velocity in an immobilised frame; n is the diffusion coefficient, representing either turbulent mixing in the water column or biodiffusion in the
Fig. 1. Coupled benthic-pelagic system sketch.
sediments. j is the porosity when Eq. (1) is applied to the dissolved tracer, and j1 in the case of sediments. To solve Eq. (1) in both the water column (zintozoztop) and sediments (zbotozozint), the location of the sediment–water interface (zint) has to be known. Assuming no fluxes through the lateral boundaries, and an immobilised, impermeable base beneath the deep sediments, the displacement of the sediment–water interface can be related to the material fluxes across the interface @zint ¼ ðF s þ F w Þ. @t
(2)
Here Fs, Fw are the normalised volumetric fluxes of solid and liquid phases across the water–sediment interface, respectively.
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
The volumetric flux of solids can be expressed as Fs ¼
Qs
rs
,
(3)
where Qs (kg/(m2 s)) is the resuspension or deposition flux of the sediments, and rs is the density of the sediment grains. To simulate liquid fluxes, we assumed that an amount of pore water expelled from sediments during resuspension is a function of the observed sediment porosity, and an amount of water trapped in sediments during deposition is a function of the porosity of fresh sediment deposits. The corresponding fluxes can be expressed as (eQs/rs), where e is a void ratio of either resuspending or settling sediments. In consolidating or swelling sediments, water exchange between benthic and pelagic layers is also influenced by an expulsion of pore water from consolidating sediments, and entrainment of water into swelling sediments. The corresponding fluxes across the sediment–water interface can as a be represented function of the interface velocity U sc z¼z , where Ucs is int the velocity of consolidating or swelling particles. Combining fluxes due to resuspension/deposition and swelling/consolidation gives a net flux of water through the sediment–water interface F w ¼ U sc z¼z
int
þ
Q s . rs
(4)
2.2. Boundary conditions Boundary flux conditions at the top of the water column and at the bottom of the sediment bed are assigned as zero for any tracer @ C ¼ 0; z ¼ ztop ; z ¼ zbot . j U˜ n (5) @z Material fluxes across the ‘‘water–sediment’’ interface are specified taking into account resuspension/deposition and swelling/consolidation @ C ¼ jq; z ¼ zint . j U˜ n (6) @z Here the flux q represents either sediment resuspension/ deposition (q ¼ Qs) or advection of the dissolved tracer across the sediment–water interface (q ¼ Fw Cd). Note that biodiffusion across water and sediments associated with bioturbation and bioirrigation could be easily incorporated into the model by adding the corresponding fluxes to q.
1169
Here Qis is the resuspension/deposition flux of the ith fraction of the particulate tracer, and ris is the density of the sediment grains, ei is the void ratio of either resuspending or settling sediments.
3. Key processes and parameters To solve (1)–(6), resuspension/deposition fluxes (Qs), settling velocity of suspended sediment (Us), velocity of consolidating particles (Ucs), vertical water velocity (Uw), and velocity of the numerical grid levels (Vg) should be specified. In this paper we utilise commonly used formulations for sediment resuspension and deposition, assuming a constant settling velocity for the suspended particles, and employing simple empirical formulations to specify sediment and pore water velocities in consolidating sediments. This section briefly outlines these formulations, and provides a more detailed description of the time-varying numerical grid specifications.
3.1. Sediment velocity There is a large body of literature concerning the settling of suspended particles in turbulent flows (e.g. Dyer, 1989; van Rijn, 1993; Thomas et al., 1999). Settling velocities are typically parameterised as a function of local hydrological characteristics and physical properties of the sediment grains. Detailed discussion of these aspects is beyond the scope of this article. For the sake of simplicity and without any loss of generality, we assume a constant settling velocity for the suspended particles. The level of accuracy and complexity in representing the velocity of bottom sediments varies from one study to another depending on data availability and modelling objectives. Geotechnical engineers have made considerable progress in understanding and predicting self-weight consolidation of fine-grained sediments (Toorman, 1996). However, modelling of natural sediments, typically represented by chemically and biologically active mixtures of water, sand, fines, and organic materials, is still in a development stage (Torfs et al., 1996; Mitchener and Torfs, 1996). In this study, to illustrate the model performance, we specify sediment velocity as a function of the sediment void ratio and assume that the void ratio in sediments gradually relaxes to a predefined ultimate void ratio
2.3. Mixture of multi-grain-size sediments
@U sc 1 @ 1 ð m Þ ¼ . ð1 þ Þ @t ð1 þ Þ T c @z
Having sediments represented by a number of size fractions (i ¼ 1n), implies that Eq. (1) with the boundary conditions (5, 6) is applied to each sediment fraction, and fluxes (3) and (4) are reformulated as follows:
Here e is the actual void ratio, em is the ultimate void ratio, and Tc is the time scale of the seabed swelling or consolidation. Note that settling velocities of consolidating particles in (7) are the same for all sediment classes. Some justification for this choice provide laboratory studies showing that in bottom sediments represented by concentrated mixtures of particles, the sediment grains tend to form a matrix structure, which is generally able to support particles without falling through (Torfs et al., 1996).
Fs ¼
n X Qs i
i¼1
rsi
(30 )
,
F w ¼ U sc z¼z
int
þ
n X i Q si . s i¼1
ri
(40 )
(7)
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
Eq. (7) represents mass conservation of a uniform medium with a time-varying porosity (Appendix A). To solve (7), the sediment velocity must be specified at the bottom of the sediment layer: Ucs ¼ 0; z ¼ zbot.
3.2. Water velocity An equation for vertical component of the water velocity in water column and in sediments is derived from the mass balance equation for the liquid phase (Appendix A) w @ jU˜ @V g @j ¼ þj . (8) @z @t @z Eq. (8) expresses water flow as a function of the timevarying porosity and sediment thickness, and thus gives an integral representation of various specific liquid fluxes, associated with sediment erosion/deposition, swelling/ consolidation, settling and bioturbation processes. To solve (8), the water velocity must be specified at the w bottom of the sediment layer: jU˜ ¼ 0; z ¼ zbot.
3.3. Resuspension and deposition fluxes Commonly used semi-empirical formulas have been employed to represent sediment resuspension and deposition on a cohesive or non-cohesive sediment bed. On a sandy seabed, it is assumed that under steady, uniform flow and sediment loading conditions, an equilibrium distribution of sediments in the water column tends to be established, with the resuspension and deposition fluxes cancelling each other. An appropriate sediment flux boundary condition can be specified using the Smith and McLean formulation (Smith and McLean, 1977). In cohesive sediments the bottom boundary condition is formulated using Ariathurai and Krone (1976), and Partheniades (1965) formulations. The sediment fluxes are functions of the bottom shear stresses, which are calculated using the Grant and Madsen approach (Madsen, 1994).
3.5. Numerical grid Time-varying numerical grids can generate apparent material fluxes between the grid layers, introducing numerical diffusion. When this artificial diffusion is comparable or exceeds diffusion due to physical or biological processes, the predicted distribution of any property in water or in sediments will be distorted. Diffusion in sediments is typically orders of magnitude less intensive than mixing in surface water. Hence, extra care has to be exercised when choosing numerical grids for sediments. The numerical grid in sediments is also supposed to have high vertical resolution near the sediment–water interface, where the sediment and water properties often exhibit high gradients. Stretched depth-adapted numerical grids are commonly used in atmospheric and ocean modelling to maintain high spatial resolution in areas of elevated variability of simulated fields (Phillips, 1957; Blumberg and Mellor, 1983). Numerical grids based on sigma coordinates, having thickness of grid layers scaled by water depth provide an example of such grid configurations. Using notations introduced in Fig. 1, a sigma grid (˜z) in sediments can be specified as follows: z˜ k ¼ za þ sk ðza zbot Þ;
sk ¼
za z˜ k za zbot
;
za 4˜zk Xzbot
(9)
za 4˜zk pzbot .
(10)
t¼t0
0.06 Height over undisturbed seabed, m
1170
Sigma grid levels 0.04 0.02 0 -0.02 0
1
2
3
4
3
4
5
Days 3.4. Active layer
Height over undisturbed seabed, m
0.06 Processes that control thickness of an active, top sediment layer (ha ¼ zintza) are not well understood. Previous methods used to prescribe the erosion depth in coastal sediments vary considerably. Reed et al. (1999) prescribed a thickness of 0.3 mm when simulating storm resuspension on the Eel River shelf (Northern California). Harris and Wiberg (2001) calculated thickness of this layer on sandy sediments as a function of mixing depth developed by migrating ripples or sheet flow. In silty sediments the active layer depth was a function of the bottom shear stress. In our paper, we assume that the thickness of the top sediment layer is constant throughout the simulation and its value has to be specified through the model calibration study.
K grid levels 0.04 0.02 0 -0.02 0
1
2
5
Days Fig. 2. Stretched numerical grids (a) sigma grid and (b) k-grid. Zero level represents initial location of sediment–water interface.
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
(11) Here D is a predefined minimum thickness of the sediment layers, and (n1) denotes one time step back in time. For the sake of convenience in this paper we will refer to grid (11) as a k-grid. During the deposition event, this grid accommodates new sediments by varying only the top sediment layers. During the resuspension event, the grid layers in deep sediments remain immobilised until the erosion depth is large enough to reach this layer. The k-grid levels in deep sediments are less mobile than sigma levels and hence less diffusive, however, as illustrated in Fig. 2b, the k-grid fails to maintain a high vertical resolution of the top sediments during the deposition event. To achieve both high vertical resolution in the top sediments and low numerical diffusion in deep layers, we combined the k-grid (11) and the sigma grid (9) into a single k-sigma grid, which is defined as follows: zk ¼ z^ k þ pð˜zk z^ k Þ;
za 4zk Xzbot .
(12)
Here p is a weight coefficient varying from zero to one. An evolution of the combined k-sigma grid levels is illustrated in Fig. 3. The sigma component of the grid adjusts the coordinate levels to long-term variations of the bathymetry, while the k-grid component enables the top sediment layer to follow high-frequency oscillations of the sediment thickness. Varying the p-value changes the time scales of the grid relaxation to the sediment thickness. Increasing p freezes the deep sediment levels, reducing numerical diffusion in deep sediments and at the same time reducing resolution of the top sediments. Decreasing the p-value improves the resolution of the top sediments, but tends to introduce higher numerical diffusion in deep layers. From (12), the velocity of the grid levels can be calculated as Vg ¼ @zk/@t. 3.6. Numerical solution For an actual simulation, the sediment concentration in water column and in sediment bed is updated, followed
Height over undisturbed seabed, m
0.05 0 -0.05 K-sigma grid levels -0.1 -0.15 0
2
4
6
8
10
Days 0 Height over undisturbed seabed, m
An evolution of the levels in this grid is illustrated in Fig. 2a. The upper curve indicates location of the top level of the sediment grid, which coincides with the sediment water interface and varies with time due to regular resuspension and deposition of sediments and due to gradual swelling of the seabed. Zero level denotes the location of the initial, undisturbed seabed surface. The sigma levels maintain refined resolution near the sediment–water interface throughout the simulation but, as shown later, they introduce large numerical diffusion in deep sediments. To reduce apparent fluxes, additional constraints must be imposed on the grid levels. The following equation provides an example of a grid, which has reduced mobility of the grid levels under the top sediments: 8 @za ðn1Þ > > o0; za 4z^ k Xzbot : < min½ðz^ kþ DÞ; z^ k ; @t z^ k ¼ @z ðn1Þ a > > X0; za 4z^ k Xzbot : : z^ k ; @t
1171
-0.03 -0.06 -0.09 -0.12 K-sigma grid levels -0.15 -0.18 0
2
4
6
8
10
Days Fig. 3. K-sigma grid levels in (a) erodible swelling and (b) erodible consolidating sediments (p ¼ 1e-2). Zero level represents initial location of sediment-water interface.
by the calculation of liquid fluxes, and solution of the transport equation for the dissolved matter. The advection-diffusion Eq. (1), is expressed in a conservative form @ @ ðC jZÞk ¼ @z j U˜ n C . (13) @t @z k Boundary conditions for (13) are provided by (5) and (6). The numerical model utilises implicit finite difference schemes to approximate advective and diffusive terms. The governing equations are transformed into a system of three-diagonal, algebraic equations and solved numerically. 4. Diffusion experiments A number of numerical experiments have been carried out to investigate numerical diffusion introduced by stretched numerical grids. Simulations involved modelling of the vertical transport of a conservative dissolved tracer in a coupled benthic-pelagic system using varying grid specifications. Resuspension and deposition were driven by tidal currents, represented by a harmonic reference velocity, which was extrapolated over the whole water depth using a log-profile approximation. Turbulent mixing in the water column was based on the level 2 turbulence closure scheme of Mellor–Yamada (1982). To test the model performance under conditions favouring the development of numerical diffusion, up to 1 cm of
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
sediment was eroded and deposited during every tidal cycle. Long-term changes in sediment thicknesses were driven by sediment swelling or consolidation processes. The dissolved tracer was initially allocated to a depth of 25 cm below the water–sediment interface (i.e. in a layer between 24 and 26 cm depth). The initial concentration of the tracer in the sediment layer was 100 mg/L and the concentration in the water column was zero. The background diffusion coefficient in sediments was constant throughout the sediment depth and equal to 0.5e9 m2/s. The initial porosity in sediments was 0.55, and final porosities were 0.65 and 0.4 for the swelling and consolidating sediments, respectively. Grid weighting parameter p (Eq. (12)) was 1e-2, unless otherwise stated. The water depth of 10 m was resolved into 10 layers. The initial thickness of sediments was 0.5 m, and sediments were resolved into 51 grid levels. The modelling period was 40 days. Simulation scenarios were subdivided into 5 groups. The first group consists of one experiment (A), which simulates molecular diffusion of the dissolved tracer in still sediments. During this simulation, sediments are neither resuspended/deposited nor consolidating/swelling, the numerical grid is immobilised and does not contribute to the numerical diffusion. All numerical grids under these conditions produce exactly the same concentration profiles. The second group of experiments deals with the diffusion in sediments that undergo tidal resuspension and deposition (B1–B3). Swelling and consolidation are switched off and there are no long-term changes in sediment thickness. The next two groups of scenarios simulate diffusion and advection in erodible sediments, with either a swelling (C1–C3) or consolidating (D1–D3) seabed. Finally, the last group of experiments (E1–E3) simulates concentration profiles in consolidating
sediments having the resuspension and deposition processes switched off. All simulations utilise the same numerical schemes, and any variations in the vertical profiles predicted by different scenarios within a single group of experiments are due to varying specifications of the numerical grids (Table 1) Fig. 4 shows vertical profiles of concentration predicted in the erodible seabed according to scenarios B1–B3. The sediments are neither swelling nor consolidating. The plot also includes the vertical profile predicted in still, unerodible sediments (scenario A), which provides a reference profile to test predictions based on time-varying grids. As shown in Fig. 4, in erodible sediments the concentrations predicted by k and k-sigma grids (scenarios B1 and B3) coincide with the concentration predicted
0 still grid (A)
Height over undisturbed seabed, m
1172
k grid (B1) sigma grid (B2)
-0.1
k-sigma grid (B3)
-0.2 -0.3 -0.4 -0.5 0
5
10 15 Concentration, mg/L
20
Fig. 4. Vertical profiles of a dissolved tracer in a still (scenario A) and erodible unconsolidating (scenarios B1–B3) sediments. Zero depth corresponds to initial location of sediment-water interface.
Table 1 Numerical experiments with a varying grid configuration. Scenario number
Resuspension/deposition
Swelling
Consolidating
Grid type
A
ks
B1 B2 B3
+ + +
s
C1 C2 C3
+ + +
+ + +
k
D1 D2 D3
+ + +
+ + +
k
E1 E2 E3
+ + +
k
Simulations involved modelling of tidally driven resuspension and deposition of sediments over consolidating/swelling seabed.
k ks
s ks
s ks
s ks
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
0.2 k grid (C1)
Height over undisturbed seabed, m
0.1
sigma grid (C2) k-sigma grid (C3)
0 -0.1 -0.2 -0.3 -0.4 -0.5
5 10 Concentration, mg/L
0
15
0 Height over undisturbed seabed, m
k grid (D1) sigma grid (D2)
-0.1
k-sigma grid (D3)
-0.2 -0.3
0 k grid (E1)
Height over undisturbed seabed, m
in still sediments (scenario A) over most of the benthic layer. A slight mismatch near the sediment surface is attributed to the varying thickness of sediments in scenario with the erodible seabed. The largest discrepancies between the reference concentrations and test profile are produced by the sigma-grid simulation (scenario B2) which grossly underestimates peak concentration in sediments. The third and the fourth groups of experiments (scenarios C1–C3 and D1–D3) simulate diffusion and advection of the dissolved tracer in erodible seabed with either swelling (Fig. 5a) or consolidating (Fig. 5b) sediments. The k-grid layers during these simulations remain immobilised below the top varying layer of sediment, and the corresponding concentrations (scenarios C1 and D1) provide reference values to test the performance of the stretched, sigma and k-sigma, grids. As shown in Fig. 5, the modelling based on k and k-sigma grids (scenarios C1–C3 and D1–D3) predicts similar concentration profiles, while sigma grids, as expected, again heavily underestimate peak concentrations (scenarios C2 and D2). Line markers on Fig. 5 indicate location of the numerical grid cell centres. In the case of a consolidating seabed (Fig. 5b), all numerical grids
1173
sigma grid (E2)
-0.1
k-sigma grid (E3)
-0.2 -0.3 -0.4 -0.5 0
5
10 15 Concentration, mg/L
20
Fig. 6. Vertical profiles of dissolved tracer in unerodible consolidating sediments.
maintain the high vertical resolution of top sediments. In swelling sediments (Fig. 5a), the k-grid resolves the top 12 cm with only 2 layers, while the k-sigma and sigma grids give an equally good resolution of the top 12 cm of sediments, resolving them into 12 layers. Resuspension and deposition of sediments enhances solute exchange between water and sediments. Fig. 6 shows vertical profiles predicted for the scenarios with consolidating sediments having sediment resuspension and deposition switched off (scenarios E1–E3). The peak of the predicted concentrations in this case reaches the surface of the sediment-water interface and the concentration levels are much higher than the values predicted by the scenario with consolidating and erodible sediments (Fig. 5b). As illustrated in Figs. 4 and 5, the quality of simulations based on sigma grids is particularly poor when sediments are subject to tidal resuspension and deposition. Having sediment resuspension and deposition switched off (Fig. 6) results in a much better agreement between the reference concentrations and concentrations predicted by the sigma-grid model (compare Figs. 6 and 5b). The low quality of the predictions in the case of erodible sediments is attributed to small-scale high-frequency oscillations of sigma levels driven by repeated resuspension and deposition of sediments. Small apparent fluxes associated with the oscillating sediment thickness when integrated over a long time introduce large numerical diffusion. Having those high-frequency oscillations filtered out leads to lower levels of the numerical diffusion and improves quality of simulations, as illustrated by simulations based on k-sigma grid.
-0.4 5. Discussion
-0.5 0
2
4 6 Concentration, mg/L
8
10
Fig. 5. Vertical profiles of dissolved tracer in (a) erodible swelling and (b) erodible consolidating sediments (scenarios C1–C3 and D1–D3, respectively).
Coupled hydrodynamic and sediment transport models often share the same numerical grid configuration in the water column. Having a numerical algorithm for the sediment transport customised to a particular grid specification (e.g. z, sigma, or isopycnal) typically implies
ARTICLE IN PRESS 1174
N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
that whenever there is the need to simulate the same processes on another grid type, the corresponding numerical algorithm has to be altered. According to the formulation developed in this study, the vertical grid structure for suspended sediments has to satisfy only a few general constraints, including alignment of the uppermost and the bottom grid levels along the surface and bottom boundaries, and monotonic continuity of the grid levels in between (the grid levels are not allowed to cross each other). Because of the high flexibility of the grid settings, the developed transport model can be coupled to hydrodynamic models having varying grid specifications without altering the numerical algorithm for the vertical transport equations. In sediments, the choice of the numerical grid is typically constrained by low background diffusion, high spatial gradients, and the dynamic nature of the sediment boundary. In this paper we provide an example of a grid that suppresses numerical diffusion in deep sediments and maintains a high vertical resolution of the top sediments by applying varying time scales to adjust the grid levels to the sediment thickness. Given a large variety of potential applications and diversity of the environmental settings, many other alternative configurations of the numerical grid settings can be designed. In some applications, one may prefer to specify grid levels in sediments as a function of a mixed depth layer, in another application it may be a function of the vertical profile with a vertical resolution dynamically improving in areas with high gradients. The grid velocity can be set equal to the velocity of consolidating sediments, thus reducing the time-varying Eulerian grid to a Lagrangian frame. The formulation developed in this paper, we believe, provides a framework for making all these changes feasible within a single numerical model. The dynamics of the numerical grid developed in this paper are a function of the weight coefficient p. Decreasing the p-value freezes the deep sediment layers, and the model adjusts to the varying geometry by changing the thickness of the top sediment layers (and by coarsening sediment resolution during deposition events). Increasing the p-value improves resolution of the top sediments, but introduces higher numerical diffusion in deep sediments. The actual choice of the p is a trade-off between the need for high resolution of the top sediments and the requirement of low numerical diffusion in deep layers. According to our observations, during long-term seasonal simulations, a p-value of 0.01 enables one to filter out highfrequency oscillations from the deep grid layers and at the same time to maintain sufficiently high resolution of the top sediments. Short-term fluctuations of the sediment thickness during such simulations are attributed to the sediment resuspension due to tidal currents, episodic high swell and storms events, etc. Long-term variations in sediment thickness are driven by the sediment swelling/ consolidation and divergence/convergence of the horizontal sediment fluxes. A first-order upwind numerical scheme was used to approximate advective terms in our simulations. Implementing higher order, less diffusive approximations would further improve the accuracy of the model
predictions. However, it is worthwhile mentioning here, that such schemes are computationally expensive. In multidisciplinary studies, where the method employed for a sub-system simulation has to be consistent with the efficiency of the whole modelling package, extra accuracy in modelling sediment profiles may not always merit losses in the simulation speed. The model developments undertaken here remain largely within the limits of the mass conservation law. Many important processes relevant to sediment dynamics, such as flocculation of fines in a turbulent flow, selfweight consolidation of cohesive sediments, development of the bed forms, bottom friction under combined waves and currents etc., were excluded from consideration. Our main objective was to develop transport equations, in a form sufficiently general for most envisaged applications of the model, which would allow a flexible choice of the sediment processes depending on the modelling objectives. Given interdependent nature of the sediment and dissolved materials transport in the coastal systems, we thought it important to incorporate into the sediment model additional equations governing transport of dissolved substances. The model described in this paper has been coupled to three-dimensional environmental modelling suite EMS (Herzfeld, 2006; Walker and Sherwood, 1997) and tested through a number of applications around the Australian coast (Robson et al., 2006; Margvelashvili et al., 2008).
6. Conclusions A one-dimensional vertical numerical model of sediment and dissolved material transport in a coupled benthic-pelagic layer has been developed. The model governing equations are formulated in a time-varying Eulerian frame and then solved numerically, which provides a consistent approach for incorporating resuspension/deposition, swelling/consolidation, and biodiffusion processes within a single modelling framework. A number of numerical experiments, with varying grid specifications, were carried out to test diffusive properties of the model. It was shown that stretched numerical grids, having varying time scales of the grid adjustment to the sediment thickness, can provide high resolution of the top sediments, while maintaining low numerical diffusion in deep layers. The developed model provides a test bed for testing scientific hypothesis concerning coastal sediment transport. It can also be used to support biogeochemical studies in a coupled benthic-pelagic system.
Acknowledgments The development of the model was supported over many years within a variety of projects including NWSJEMS, Ord Bonaparte Program, SRFME, CRC for Coastal Zone, Estuary & Waterway Management, CRC Reef Torres Strait program and Aquafin CRC.
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
Appendix A. Mass balance equations In an Eulerian immobilised coordinate frame z, a general mass balance equation for the particulate or dissolved tracer can be expressed as follows (Boudreau, 1997): @jC @ðjCUÞ @ @C þ ¼ . (A.1) jn @t @z @z @z Here C is either the volumetric concentration of the particulate tracer or concentration of the dissolved tracer in pore water, U is the material velocity given in a timeinvariant coordinate frame, and j is the porosity when C is a dissolved tracer and j1 when (1) is applied to a particulate fraction. Introducing a new set of independent variables t ¼ t; z ¼ z ðtÞ;
@z ¼1 @z
(A.2)
and applying a chain rule to link derivatives in the old system to those in the new system @jC @jC @jC @z @jC @jC ¼ ¼ ; þ @t @t @z @t @z @z the mass balance Eq. (A.1) is transformed into @jC @jC @z @jCU @ @C þ ¼ jn . þ @t @z @t @z @z @z
(A.3)
(A.5)
(A.6)
and suppressing the star superscripts in the notation of the independent variables, equation for the mass conservation in a time-varying coordinate frame can be expressed as @jC @ @ @V g ¼ C jC . (A.7) j U˜ n @t @z @z @z Here U˜ ¼ U V g is the material velocity expressed in a time-varying frame. Assuming C ¼ 1, and zero diffusion coefficient in (A.7) gives a mass conservation equation for the liquid phase w
@j @ðjU˜ Þ @V g þj þ ¼ 0. @z @t @z Equation for the concentration of solids reads @V g @C @ @ ¼ U˜ n C C . @t @z @z @z
Introducing a time-varying coordinate frame Vg ¼ U, and substituting in (A.10) porosity by void ratio (1j) ¼ (1/(1+e)), reduces (A.10) to @ @U . ¼ ð1 þ Þ @z @t
(A.11)
Using the expression for the sediment velocity in a cohesive consolidating sediment (Toorman, 1996) K @se @ K h ðrs rw Þ U ¼ hw (A.12) gr @ @z ð1 þ Þ rw and introducing new coordinates @x/@z ¼ (1/(1+e)), transforms (A.11) to the finite strain consolidation equation derived by Gibson et al. (1967) @ @ Kh @se @ K h ðrs rw Þ ¼ . w w @t @x ð1 þ Þg r @ @x ð1 þ Þ r (A.13)
(A.4)
Introducing the velocity of new coordinate levels as seen by an observer in an immobilised coordinate system (note that @z*/@t is velocity of the old frame points with regard to new coordinate levels) @z , @t
side of (A.9) by the sediment density, and neglecting diffusion, transforms (A.9) into
@V g @ð1 jÞ @ ð1 jÞU˜ þ þ ð1 jÞ ¼ 0. (A.10) @t @z @z
Here se is an effective stress, and Kh is hydraulic conductivity.
After a number of manipulations involving expansion of partial derivatives and combining different terms, (A.4) can be transformed to @jC @ @z @ @ @z þ U n C jC ¼ 0. j þ @t @z @z @z @t @t
Vg ¼
1175
(A.8)
(A.9)
Having sediment represented by single class of particles: (C/r) ¼ (1j), and dividing the left and right hand
References Ariathurai, R., Krone, R.B., 1976. Finite element model for cohesive sediment transport. Journal of the Hydraulics Division 102, 323–338. Blumberg, A., Mellor, G., 1983. Diagnostic and prognostic numerical circulation studies of the South Atlantic Bight. Journal of Geophysical Research 88, 4579–4592. Boudreau, B.P., 1997. Diagenetic Models and Their Implementation: Modelling Transport and Reactions in Aquatic Sediments. Springer, Berlin, 414pp. Dyer, K.R., 1989. Sediment processes in estuaries: future research requirements. Journal of Geophysical Research 94 (C10), 14327–14339. Gibson, R.E., England, G.L., Hussey, M.L., 1967. The theory of onedimensional consolidation of saturated clays-I. Geotechnique 17, 261–273. Govindaraju, R., Ramireddygari, R., Shrestha, P., Roig, L., 1999. Continuum bed model for estuarine sediments based on nonlinear consolidation theory. Journal of Hydraulic Engineering 125, 300–304. Harris, C., Wiberg, P., 2001. A two-dimensional, time-dependent model of suspended sediment transport and bed reworking for continental shelves. Computers & Geosciences 27, 675–690. Herzfeld, M., 2006. An alternative coordinate system for solving finite difference ocean models. Ocean Modelling 14, 174–196. Luff, R., Moll, A., 2004. Seasonal dynamics of the North Sea sediments using a three-dimensional coupled sediment–water model system. Continental Shelf Research 24, 1099–1127. Madsen, O.S., 1994. Spectral wave-current bottom boundary layer flows. In: Proceedings 24 International Conference on Coastal Engineering. American Society of Civil Engineers, Kobe, Japan, pp. 384–398. Margvelashvili, N., Saint-Cast, F., Condie, S., 2008. Numerical modelling of the suspended sediment transport in Torres Strait. Continental Shelf Research 28, 2241–2256. Mellor, G., Yamada, T., 1982. Development of a turbulence closure model for geophysical fluid problems. Reviews of Geophysics and Space Physics 20, 851–875. Mitchener, H., Torfs, H., 1996. Erosion of mud/sand mixtures. Coastal Engineering 29, 1–25. Narasimhan, M., 1992. Principles of Continuum Mechanics. Wiley, New York, 584pp. Partheniades, E., 1965. Erosion and deposition of cohesive soils. Journal of the Hydraulic Division 91, 105–139. Phillips, N.A., 1957. A coordinate system having some special advantages for numerical forecasting. Journal of Meteorology 14, 184–185.
ARTICLE IN PRESS 1176
N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
Reed, C.W., Neidiroda, A.W., Swift, J.P., 1999. Modelling sediment entrainment and transport processes limited by bed armoring. Marine Geology 154, 143–154. Robson, B., Webster, I., Margvelashvili, N., Herzfeld, M., 2006. Scenario modelling: simulating the downstream effects of changes in catchment land use. Cooperative Research Centre for Coastal Zone, Estuary and Water Way Management, Technical Report No.41, Indooroopilly, Australia, 34pp. Skogen, M., Svendsen, E., Bernsten, J., Aksnes, D., Ulvestad, K., 1995. Modelling the primary production in the North Sea using a coupled 3-dimensional physical-chemical-biological ocean model. Estuarine, Coastal and Shelf Sciences 41, 545–565. Smith, J.D., McLean, S.R., 1977. Spatially averaged flow over a wavy bed. Journal of Geophysical Research 82, 1735–1746.
Thomas, D.N., Judd, S.J., Fawcett, N., 1999. Flocculation modelling: a review. Water Resources 33, 1579–1592. Toorman, E., 1996. Sedimentation and self weight consolidation: general unifying theory. Geotechnique 46, 103–113. Torfs, H., Mitchener, H., Huysentruy, H., Toorman, E., 1996. Settling and consolidation of mud/sand mixtures. Coastal Engineering 29, 27–45. Wainright, S., Hopkinson, C., 1997. Effects of sediment resuspension on organic matter processing in coastal environments: a simulation model. Journal of Marine Systems 11, 353–368. Walker, S.J., Sherwood, C.R., 1997. A transport model of Port Phillip Bay. CSIRO Division of Oceanography, Technical Report No. 39, Hobart, Australia, 59pp. Van Rijn, L., 1993. Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas. Aqua Publications, Amsterdam (Multiple pagination).
Contents lists available at ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
Stretched Eulerian coordinate model of coastal sediment transport N. Margvelashvili CSIRO Marine and Atmospheric Research, GPO Box 1538, Hobart, Tasmania 7001, Australia
a r t i c l e in fo
abstract
Article history: Received 16 January 2008 Accepted 28 March 2008
A one-dimensional vertical model of sediment and dissolved material transport in a coupled, vertically resolved benthic-pelagic system has been developed. An apparently unique feature of the model is that the governing equations are first formulated in a time-varying Eulerian frame and then solved numerically, which provides a consistent approach for incorporating resuspension/deposition, swelling/consolidation, and biodiffusion processes within a single modelling framework. Numerical experiments were carried out with varying grids to test diffusive properties of the model. A stretched numerical grid that maintains high vertical resolution of the top sediments and constrains numerical diffusion in deep layers is an integral part of the model. The developed model provides a test bed for testing scientific hypothesis pertaining to coastal sediment transport. The model can also be used to support biogeochemical studies in a coupled benthic-pelagic system. Crown Copyright & 2008 Published by Elsevier Ltd. All rights reserved.
Keywords: Coastal Sediment Numerical Model Benthic Pelagic
1. Introduction Vertical profiles of physical, chemical and biological properties in coastal water and sediments often exhibit strong spatial gradients near the sediment–water interface (Dyer, 1989; Boudreau, 1997). The location of the interface varies with time due to resuspension and deposition or swelling/consolidation of the seabed, and the displacement of this boundary may exceed typical scales of spatial variability of the vertical profiles. The dynamic nature and high spatial variability of the nearbottom fields provide a challenge for adequate description of benthic fluxes across the water and sediment layers, which is essential for better understanding and efficient management of coastal ecosystems. To address this challenge, we developed a modelling framework for coastal material transport that accommodates key sediment processes and reconciles strong vertical gradients with the dynamic nature of coastal sediments.
Tel.: +61 3 62325142; fax: +61 3 62325123.
E-mail address: [email protected]
Integrated coastal models typically include ecological reaction models and physical transport models, operating in a coupled benthic-pelagic system (Skogen et al., 1995; Wainright and Hopkinson, 1997; Luff and Moll, 2004). Transport models in sediments typically employ either a time-invariant Eulerian frame or Lagrangian coordinates. Lagrangian coordinates can provide high vertical resolution of the uppermost sediments during the simulation (Gibson et al., 1967; Govindaraju et al., 1999). However these coordinates are not well suited for simulating mixtures of sand, fines and dissolved materials (which is typically the case with natural sediments), since in a Lagrangian frame the coordinates are fixed to a given material parcel which should always embrace the same solid or dissolved particles, a particular choice being dependent on the simulation objectives (Narasimhan, 1992). In Eulerian coordinates, having all grid layers fixed would imply constant thickness of sediments. To adjust a numerical grid to the time-varying thickness, in Eulerian coordinates either new sediment layers must be added on top of existing layers or old layers must be diluted during the course of simulation. To maintain a high speed in the simulation, the resolution of top sediments must vary
0098-3004/$ - see front matter Crown Copyright & 2008 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2008.03.016
ARTICLE IN PRESS 1168
N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
dramatically with time, compromising the accuracy of the calculations. An alternative approach for simulating material transport in sediments is to solve transport equations using Eulerian grids, where all grid layers vary with time (Harris and Wiberg, 2001; Reed et al., 1999). Since the time-varying grid layers introduce large numerical diffusion in sediments, these models simulate suspended sediment transport over a two-layered seabed assuming a fully mixed environment under the top sediments. The model described in this paper extends the latter formulation to multilayered sediments. It employs timevarying Eulerian grids plus a user-specified dynamics for the grid layers. To maintain high resolution in the top sediments without sacrificing computational efficiency, the model utilises varying time scales of the grid adjustment to the sediment thickness. The developed model conserves masses of solid, liquid and dissolved phases, and provides a consistent approach for incorporating resuspension/deposition, swelling/consolidation, and biodiffusion processes within a single modelling framework.
2. Transport equations Here and throughout the paper, solid and liquid phases are assumed to be incompressible. The vertical axis of the coordinate frame is directed upward. Sediments and dissolved tracers are represented by vertical profiles of the sediment concentration Cs(z, t) (kg/(m3 mixture)), and the concentration of the dissolved tracer Cd(z, t) (kg/(m3 water)). The superscripts (s), (w), and (d) denote the solid phase, liquid phase (water), and dissolved phase, respectively. The subscripts (i) and (j) are introduced whenever sediments or dissolved traces are represented by a number of classes. All classes of particles in sediments have the same vertical velocity associated with the sediment consolidation UisUcs(z, t) and the vertical velocity of all dissolved tracers is equal to the water velocity UidUw(z, t). Sediments and water column are represented by discrete, time-varying layers of thickness Zk ¼ z(k+)z(k), with the subscript k indicating the layer number, and (k+) and (k) indexes indicating the top and bottom interfaces between layers, respectively (Fig. 1).
2.1. Mass conservation A general mass balance equation for sediments (C ¼ Cs) or dissolved tracer (C ¼ Cd), in a stretched Eulerian frame z ¼ z(t), can be expressed as follows (Appendix A): @V g @ @ @ ðC jÞ ¼ C jC . (1) j U˜ n @t @z @z @z Here Vg is the velocity of the time-varying coordinate levels; U˜ ¼ ðU V g Þ is the velocity of either a particulate or dissolved tracer expressed in a time-varying frame z ¼ z(t); U is the tracer velocity in an immobilised frame; n is the diffusion coefficient, representing either turbulent mixing in the water column or biodiffusion in the
Fig. 1. Coupled benthic-pelagic system sketch.
sediments. j is the porosity when Eq. (1) is applied to the dissolved tracer, and j1 in the case of sediments. To solve Eq. (1) in both the water column (zintozoztop) and sediments (zbotozozint), the location of the sediment–water interface (zint) has to be known. Assuming no fluxes through the lateral boundaries, and an immobilised, impermeable base beneath the deep sediments, the displacement of the sediment–water interface can be related to the material fluxes across the interface @zint ¼ ðF s þ F w Þ. @t
(2)
Here Fs, Fw are the normalised volumetric fluxes of solid and liquid phases across the water–sediment interface, respectively.
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
The volumetric flux of solids can be expressed as Fs ¼
Qs
rs
,
(3)
where Qs (kg/(m2 s)) is the resuspension or deposition flux of the sediments, and rs is the density of the sediment grains. To simulate liquid fluxes, we assumed that an amount of pore water expelled from sediments during resuspension is a function of the observed sediment porosity, and an amount of water trapped in sediments during deposition is a function of the porosity of fresh sediment deposits. The corresponding fluxes can be expressed as (eQs/rs), where e is a void ratio of either resuspending or settling sediments. In consolidating or swelling sediments, water exchange between benthic and pelagic layers is also influenced by an expulsion of pore water from consolidating sediments, and entrainment of water into swelling sediments. The corresponding fluxes across the sediment–water interface can as a be represented function of the interface velocity U sc z¼z , where Ucs is int the velocity of consolidating or swelling particles. Combining fluxes due to resuspension/deposition and swelling/consolidation gives a net flux of water through the sediment–water interface F w ¼ U sc z¼z
int
þ
Q s . rs
(4)
2.2. Boundary conditions Boundary flux conditions at the top of the water column and at the bottom of the sediment bed are assigned as zero for any tracer @ C ¼ 0; z ¼ ztop ; z ¼ zbot . j U˜ n (5) @z Material fluxes across the ‘‘water–sediment’’ interface are specified taking into account resuspension/deposition and swelling/consolidation @ C ¼ jq; z ¼ zint . j U˜ n (6) @z Here the flux q represents either sediment resuspension/ deposition (q ¼ Qs) or advection of the dissolved tracer across the sediment–water interface (q ¼ Fw Cd). Note that biodiffusion across water and sediments associated with bioturbation and bioirrigation could be easily incorporated into the model by adding the corresponding fluxes to q.
1169
Here Qis is the resuspension/deposition flux of the ith fraction of the particulate tracer, and ris is the density of the sediment grains, ei is the void ratio of either resuspending or settling sediments.
3. Key processes and parameters To solve (1)–(6), resuspension/deposition fluxes (Qs), settling velocity of suspended sediment (Us), velocity of consolidating particles (Ucs), vertical water velocity (Uw), and velocity of the numerical grid levels (Vg) should be specified. In this paper we utilise commonly used formulations for sediment resuspension and deposition, assuming a constant settling velocity for the suspended particles, and employing simple empirical formulations to specify sediment and pore water velocities in consolidating sediments. This section briefly outlines these formulations, and provides a more detailed description of the time-varying numerical grid specifications.
3.1. Sediment velocity There is a large body of literature concerning the settling of suspended particles in turbulent flows (e.g. Dyer, 1989; van Rijn, 1993; Thomas et al., 1999). Settling velocities are typically parameterised as a function of local hydrological characteristics and physical properties of the sediment grains. Detailed discussion of these aspects is beyond the scope of this article. For the sake of simplicity and without any loss of generality, we assume a constant settling velocity for the suspended particles. The level of accuracy and complexity in representing the velocity of bottom sediments varies from one study to another depending on data availability and modelling objectives. Geotechnical engineers have made considerable progress in understanding and predicting self-weight consolidation of fine-grained sediments (Toorman, 1996). However, modelling of natural sediments, typically represented by chemically and biologically active mixtures of water, sand, fines, and organic materials, is still in a development stage (Torfs et al., 1996; Mitchener and Torfs, 1996). In this study, to illustrate the model performance, we specify sediment velocity as a function of the sediment void ratio and assume that the void ratio in sediments gradually relaxes to a predefined ultimate void ratio
2.3. Mixture of multi-grain-size sediments
@U sc 1 @ 1 ð m Þ ¼ . ð1 þ Þ @t ð1 þ Þ T c @z
Having sediments represented by a number of size fractions (i ¼ 1n), implies that Eq. (1) with the boundary conditions (5, 6) is applied to each sediment fraction, and fluxes (3) and (4) are reformulated as follows:
Here e is the actual void ratio, em is the ultimate void ratio, and Tc is the time scale of the seabed swelling or consolidation. Note that settling velocities of consolidating particles in (7) are the same for all sediment classes. Some justification for this choice provide laboratory studies showing that in bottom sediments represented by concentrated mixtures of particles, the sediment grains tend to form a matrix structure, which is generally able to support particles without falling through (Torfs et al., 1996).
Fs ¼
n X Qs i
i¼1
rsi
(30 )
,
F w ¼ U sc z¼z
int
þ
n X i Q si . s i¼1
ri
(40 )
(7)
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
Eq. (7) represents mass conservation of a uniform medium with a time-varying porosity (Appendix A). To solve (7), the sediment velocity must be specified at the bottom of the sediment layer: Ucs ¼ 0; z ¼ zbot.
3.2. Water velocity An equation for vertical component of the water velocity in water column and in sediments is derived from the mass balance equation for the liquid phase (Appendix A) w @ jU˜ @V g @j ¼ þj . (8) @z @t @z Eq. (8) expresses water flow as a function of the timevarying porosity and sediment thickness, and thus gives an integral representation of various specific liquid fluxes, associated with sediment erosion/deposition, swelling/ consolidation, settling and bioturbation processes. To solve (8), the water velocity must be specified at the w bottom of the sediment layer: jU˜ ¼ 0; z ¼ zbot.
3.3. Resuspension and deposition fluxes Commonly used semi-empirical formulas have been employed to represent sediment resuspension and deposition on a cohesive or non-cohesive sediment bed. On a sandy seabed, it is assumed that under steady, uniform flow and sediment loading conditions, an equilibrium distribution of sediments in the water column tends to be established, with the resuspension and deposition fluxes cancelling each other. An appropriate sediment flux boundary condition can be specified using the Smith and McLean formulation (Smith and McLean, 1977). In cohesive sediments the bottom boundary condition is formulated using Ariathurai and Krone (1976), and Partheniades (1965) formulations. The sediment fluxes are functions of the bottom shear stresses, which are calculated using the Grant and Madsen approach (Madsen, 1994).
3.5. Numerical grid Time-varying numerical grids can generate apparent material fluxes between the grid layers, introducing numerical diffusion. When this artificial diffusion is comparable or exceeds diffusion due to physical or biological processes, the predicted distribution of any property in water or in sediments will be distorted. Diffusion in sediments is typically orders of magnitude less intensive than mixing in surface water. Hence, extra care has to be exercised when choosing numerical grids for sediments. The numerical grid in sediments is also supposed to have high vertical resolution near the sediment–water interface, where the sediment and water properties often exhibit high gradients. Stretched depth-adapted numerical grids are commonly used in atmospheric and ocean modelling to maintain high spatial resolution in areas of elevated variability of simulated fields (Phillips, 1957; Blumberg and Mellor, 1983). Numerical grids based on sigma coordinates, having thickness of grid layers scaled by water depth provide an example of such grid configurations. Using notations introduced in Fig. 1, a sigma grid (˜z) in sediments can be specified as follows: z˜ k ¼ za þ sk ðza zbot Þ;
sk ¼
za z˜ k za zbot
;
za 4˜zk Xzbot
(9)
za 4˜zk pzbot .
(10)
t¼t0
0.06 Height over undisturbed seabed, m
1170
Sigma grid levels 0.04 0.02 0 -0.02 0
1
2
3
4
3
4
5
Days 3.4. Active layer
Height over undisturbed seabed, m
0.06 Processes that control thickness of an active, top sediment layer (ha ¼ zintza) are not well understood. Previous methods used to prescribe the erosion depth in coastal sediments vary considerably. Reed et al. (1999) prescribed a thickness of 0.3 mm when simulating storm resuspension on the Eel River shelf (Northern California). Harris and Wiberg (2001) calculated thickness of this layer on sandy sediments as a function of mixing depth developed by migrating ripples or sheet flow. In silty sediments the active layer depth was a function of the bottom shear stress. In our paper, we assume that the thickness of the top sediment layer is constant throughout the simulation and its value has to be specified through the model calibration study.
K grid levels 0.04 0.02 0 -0.02 0
1
2
5
Days Fig. 2. Stretched numerical grids (a) sigma grid and (b) k-grid. Zero level represents initial location of sediment–water interface.
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
(11) Here D is a predefined minimum thickness of the sediment layers, and (n1) denotes one time step back in time. For the sake of convenience in this paper we will refer to grid (11) as a k-grid. During the deposition event, this grid accommodates new sediments by varying only the top sediment layers. During the resuspension event, the grid layers in deep sediments remain immobilised until the erosion depth is large enough to reach this layer. The k-grid levels in deep sediments are less mobile than sigma levels and hence less diffusive, however, as illustrated in Fig. 2b, the k-grid fails to maintain a high vertical resolution of the top sediments during the deposition event. To achieve both high vertical resolution in the top sediments and low numerical diffusion in deep layers, we combined the k-grid (11) and the sigma grid (9) into a single k-sigma grid, which is defined as follows: zk ¼ z^ k þ pð˜zk z^ k Þ;
za 4zk Xzbot .
(12)
Here p is a weight coefficient varying from zero to one. An evolution of the combined k-sigma grid levels is illustrated in Fig. 3. The sigma component of the grid adjusts the coordinate levels to long-term variations of the bathymetry, while the k-grid component enables the top sediment layer to follow high-frequency oscillations of the sediment thickness. Varying the p-value changes the time scales of the grid relaxation to the sediment thickness. Increasing p freezes the deep sediment levels, reducing numerical diffusion in deep sediments and at the same time reducing resolution of the top sediments. Decreasing the p-value improves the resolution of the top sediments, but tends to introduce higher numerical diffusion in deep layers. From (12), the velocity of the grid levels can be calculated as Vg ¼ @zk/@t. 3.6. Numerical solution For an actual simulation, the sediment concentration in water column and in sediment bed is updated, followed
Height over undisturbed seabed, m
0.05 0 -0.05 K-sigma grid levels -0.1 -0.15 0
2
4
6
8
10
Days 0 Height over undisturbed seabed, m
An evolution of the levels in this grid is illustrated in Fig. 2a. The upper curve indicates location of the top level of the sediment grid, which coincides with the sediment water interface and varies with time due to regular resuspension and deposition of sediments and due to gradual swelling of the seabed. Zero level denotes the location of the initial, undisturbed seabed surface. The sigma levels maintain refined resolution near the sediment–water interface throughout the simulation but, as shown later, they introduce large numerical diffusion in deep sediments. To reduce apparent fluxes, additional constraints must be imposed on the grid levels. The following equation provides an example of a grid, which has reduced mobility of the grid levels under the top sediments: 8 @za ðn1Þ > > o0; za 4z^ k Xzbot : < min½ðz^ kþ DÞ; z^ k ; @t z^ k ¼ @z ðn1Þ a > > X0; za 4z^ k Xzbot : : z^ k ; @t
1171
-0.03 -0.06 -0.09 -0.12 K-sigma grid levels -0.15 -0.18 0
2
4
6
8
10
Days Fig. 3. K-sigma grid levels in (a) erodible swelling and (b) erodible consolidating sediments (p ¼ 1e-2). Zero level represents initial location of sediment-water interface.
by the calculation of liquid fluxes, and solution of the transport equation for the dissolved matter. The advection-diffusion Eq. (1), is expressed in a conservative form @ @ ðC jZÞk ¼ @z j U˜ n C . (13) @t @z k Boundary conditions for (13) are provided by (5) and (6). The numerical model utilises implicit finite difference schemes to approximate advective and diffusive terms. The governing equations are transformed into a system of three-diagonal, algebraic equations and solved numerically. 4. Diffusion experiments A number of numerical experiments have been carried out to investigate numerical diffusion introduced by stretched numerical grids. Simulations involved modelling of the vertical transport of a conservative dissolved tracer in a coupled benthic-pelagic system using varying grid specifications. Resuspension and deposition were driven by tidal currents, represented by a harmonic reference velocity, which was extrapolated over the whole water depth using a log-profile approximation. Turbulent mixing in the water column was based on the level 2 turbulence closure scheme of Mellor–Yamada (1982). To test the model performance under conditions favouring the development of numerical diffusion, up to 1 cm of
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
sediment was eroded and deposited during every tidal cycle. Long-term changes in sediment thicknesses were driven by sediment swelling or consolidation processes. The dissolved tracer was initially allocated to a depth of 25 cm below the water–sediment interface (i.e. in a layer between 24 and 26 cm depth). The initial concentration of the tracer in the sediment layer was 100 mg/L and the concentration in the water column was zero. The background diffusion coefficient in sediments was constant throughout the sediment depth and equal to 0.5e9 m2/s. The initial porosity in sediments was 0.55, and final porosities were 0.65 and 0.4 for the swelling and consolidating sediments, respectively. Grid weighting parameter p (Eq. (12)) was 1e-2, unless otherwise stated. The water depth of 10 m was resolved into 10 layers. The initial thickness of sediments was 0.5 m, and sediments were resolved into 51 grid levels. The modelling period was 40 days. Simulation scenarios were subdivided into 5 groups. The first group consists of one experiment (A), which simulates molecular diffusion of the dissolved tracer in still sediments. During this simulation, sediments are neither resuspended/deposited nor consolidating/swelling, the numerical grid is immobilised and does not contribute to the numerical diffusion. All numerical grids under these conditions produce exactly the same concentration profiles. The second group of experiments deals with the diffusion in sediments that undergo tidal resuspension and deposition (B1–B3). Swelling and consolidation are switched off and there are no long-term changes in sediment thickness. The next two groups of scenarios simulate diffusion and advection in erodible sediments, with either a swelling (C1–C3) or consolidating (D1–D3) seabed. Finally, the last group of experiments (E1–E3) simulates concentration profiles in consolidating
sediments having the resuspension and deposition processes switched off. All simulations utilise the same numerical schemes, and any variations in the vertical profiles predicted by different scenarios within a single group of experiments are due to varying specifications of the numerical grids (Table 1) Fig. 4 shows vertical profiles of concentration predicted in the erodible seabed according to scenarios B1–B3. The sediments are neither swelling nor consolidating. The plot also includes the vertical profile predicted in still, unerodible sediments (scenario A), which provides a reference profile to test predictions based on time-varying grids. As shown in Fig. 4, in erodible sediments the concentrations predicted by k and k-sigma grids (scenarios B1 and B3) coincide with the concentration predicted
0 still grid (A)
Height over undisturbed seabed, m
1172
k grid (B1) sigma grid (B2)
-0.1
k-sigma grid (B3)
-0.2 -0.3 -0.4 -0.5 0
5
10 15 Concentration, mg/L
20
Fig. 4. Vertical profiles of a dissolved tracer in a still (scenario A) and erodible unconsolidating (scenarios B1–B3) sediments. Zero depth corresponds to initial location of sediment-water interface.
Table 1 Numerical experiments with a varying grid configuration. Scenario number
Resuspension/deposition
Swelling
Consolidating
Grid type
A
ks
B1 B2 B3
+ + +
s
C1 C2 C3
+ + +
+ + +
k
D1 D2 D3
+ + +
+ + +
k
E1 E2 E3
+ + +
k
Simulations involved modelling of tidally driven resuspension and deposition of sediments over consolidating/swelling seabed.
k ks
s ks
s ks
s ks
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
0.2 k grid (C1)
Height over undisturbed seabed, m
0.1
sigma grid (C2) k-sigma grid (C3)
0 -0.1 -0.2 -0.3 -0.4 -0.5
5 10 Concentration, mg/L
0
15
0 Height over undisturbed seabed, m
k grid (D1) sigma grid (D2)
-0.1
k-sigma grid (D3)
-0.2 -0.3
0 k grid (E1)
Height over undisturbed seabed, m
in still sediments (scenario A) over most of the benthic layer. A slight mismatch near the sediment surface is attributed to the varying thickness of sediments in scenario with the erodible seabed. The largest discrepancies between the reference concentrations and test profile are produced by the sigma-grid simulation (scenario B2) which grossly underestimates peak concentration in sediments. The third and the fourth groups of experiments (scenarios C1–C3 and D1–D3) simulate diffusion and advection of the dissolved tracer in erodible seabed with either swelling (Fig. 5a) or consolidating (Fig. 5b) sediments. The k-grid layers during these simulations remain immobilised below the top varying layer of sediment, and the corresponding concentrations (scenarios C1 and D1) provide reference values to test the performance of the stretched, sigma and k-sigma, grids. As shown in Fig. 5, the modelling based on k and k-sigma grids (scenarios C1–C3 and D1–D3) predicts similar concentration profiles, while sigma grids, as expected, again heavily underestimate peak concentrations (scenarios C2 and D2). Line markers on Fig. 5 indicate location of the numerical grid cell centres. In the case of a consolidating seabed (Fig. 5b), all numerical grids
1173
sigma grid (E2)
-0.1
k-sigma grid (E3)
-0.2 -0.3 -0.4 -0.5 0
5
10 15 Concentration, mg/L
20
Fig. 6. Vertical profiles of dissolved tracer in unerodible consolidating sediments.
maintain the high vertical resolution of top sediments. In swelling sediments (Fig. 5a), the k-grid resolves the top 12 cm with only 2 layers, while the k-sigma and sigma grids give an equally good resolution of the top 12 cm of sediments, resolving them into 12 layers. Resuspension and deposition of sediments enhances solute exchange between water and sediments. Fig. 6 shows vertical profiles predicted for the scenarios with consolidating sediments having sediment resuspension and deposition switched off (scenarios E1–E3). The peak of the predicted concentrations in this case reaches the surface of the sediment-water interface and the concentration levels are much higher than the values predicted by the scenario with consolidating and erodible sediments (Fig. 5b). As illustrated in Figs. 4 and 5, the quality of simulations based on sigma grids is particularly poor when sediments are subject to tidal resuspension and deposition. Having sediment resuspension and deposition switched off (Fig. 6) results in a much better agreement between the reference concentrations and concentrations predicted by the sigma-grid model (compare Figs. 6 and 5b). The low quality of the predictions in the case of erodible sediments is attributed to small-scale high-frequency oscillations of sigma levels driven by repeated resuspension and deposition of sediments. Small apparent fluxes associated with the oscillating sediment thickness when integrated over a long time introduce large numerical diffusion. Having those high-frequency oscillations filtered out leads to lower levels of the numerical diffusion and improves quality of simulations, as illustrated by simulations based on k-sigma grid.
-0.4 5. Discussion
-0.5 0
2
4 6 Concentration, mg/L
8
10
Fig. 5. Vertical profiles of dissolved tracer in (a) erodible swelling and (b) erodible consolidating sediments (scenarios C1–C3 and D1–D3, respectively).
Coupled hydrodynamic and sediment transport models often share the same numerical grid configuration in the water column. Having a numerical algorithm for the sediment transport customised to a particular grid specification (e.g. z, sigma, or isopycnal) typically implies
ARTICLE IN PRESS 1174
N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
that whenever there is the need to simulate the same processes on another grid type, the corresponding numerical algorithm has to be altered. According to the formulation developed in this study, the vertical grid structure for suspended sediments has to satisfy only a few general constraints, including alignment of the uppermost and the bottom grid levels along the surface and bottom boundaries, and monotonic continuity of the grid levels in between (the grid levels are not allowed to cross each other). Because of the high flexibility of the grid settings, the developed transport model can be coupled to hydrodynamic models having varying grid specifications without altering the numerical algorithm for the vertical transport equations. In sediments, the choice of the numerical grid is typically constrained by low background diffusion, high spatial gradients, and the dynamic nature of the sediment boundary. In this paper we provide an example of a grid that suppresses numerical diffusion in deep sediments and maintains a high vertical resolution of the top sediments by applying varying time scales to adjust the grid levels to the sediment thickness. Given a large variety of potential applications and diversity of the environmental settings, many other alternative configurations of the numerical grid settings can be designed. In some applications, one may prefer to specify grid levels in sediments as a function of a mixed depth layer, in another application it may be a function of the vertical profile with a vertical resolution dynamically improving in areas with high gradients. The grid velocity can be set equal to the velocity of consolidating sediments, thus reducing the time-varying Eulerian grid to a Lagrangian frame. The formulation developed in this paper, we believe, provides a framework for making all these changes feasible within a single numerical model. The dynamics of the numerical grid developed in this paper are a function of the weight coefficient p. Decreasing the p-value freezes the deep sediment layers, and the model adjusts to the varying geometry by changing the thickness of the top sediment layers (and by coarsening sediment resolution during deposition events). Increasing the p-value improves resolution of the top sediments, but introduces higher numerical diffusion in deep sediments. The actual choice of the p is a trade-off between the need for high resolution of the top sediments and the requirement of low numerical diffusion in deep layers. According to our observations, during long-term seasonal simulations, a p-value of 0.01 enables one to filter out highfrequency oscillations from the deep grid layers and at the same time to maintain sufficiently high resolution of the top sediments. Short-term fluctuations of the sediment thickness during such simulations are attributed to the sediment resuspension due to tidal currents, episodic high swell and storms events, etc. Long-term variations in sediment thickness are driven by the sediment swelling/ consolidation and divergence/convergence of the horizontal sediment fluxes. A first-order upwind numerical scheme was used to approximate advective terms in our simulations. Implementing higher order, less diffusive approximations would further improve the accuracy of the model
predictions. However, it is worthwhile mentioning here, that such schemes are computationally expensive. In multidisciplinary studies, where the method employed for a sub-system simulation has to be consistent with the efficiency of the whole modelling package, extra accuracy in modelling sediment profiles may not always merit losses in the simulation speed. The model developments undertaken here remain largely within the limits of the mass conservation law. Many important processes relevant to sediment dynamics, such as flocculation of fines in a turbulent flow, selfweight consolidation of cohesive sediments, development of the bed forms, bottom friction under combined waves and currents etc., were excluded from consideration. Our main objective was to develop transport equations, in a form sufficiently general for most envisaged applications of the model, which would allow a flexible choice of the sediment processes depending on the modelling objectives. Given interdependent nature of the sediment and dissolved materials transport in the coastal systems, we thought it important to incorporate into the sediment model additional equations governing transport of dissolved substances. The model described in this paper has been coupled to three-dimensional environmental modelling suite EMS (Herzfeld, 2006; Walker and Sherwood, 1997) and tested through a number of applications around the Australian coast (Robson et al., 2006; Margvelashvili et al., 2008).
6. Conclusions A one-dimensional vertical numerical model of sediment and dissolved material transport in a coupled benthic-pelagic layer has been developed. The model governing equations are formulated in a time-varying Eulerian frame and then solved numerically, which provides a consistent approach for incorporating resuspension/deposition, swelling/consolidation, and biodiffusion processes within a single modelling framework. A number of numerical experiments, with varying grid specifications, were carried out to test diffusive properties of the model. It was shown that stretched numerical grids, having varying time scales of the grid adjustment to the sediment thickness, can provide high resolution of the top sediments, while maintaining low numerical diffusion in deep layers. The developed model provides a test bed for testing scientific hypothesis concerning coastal sediment transport. It can also be used to support biogeochemical studies in a coupled benthic-pelagic system.
Acknowledgments The development of the model was supported over many years within a variety of projects including NWSJEMS, Ord Bonaparte Program, SRFME, CRC for Coastal Zone, Estuary & Waterway Management, CRC Reef Torres Strait program and Aquafin CRC.
ARTICLE IN PRESS N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
Appendix A. Mass balance equations In an Eulerian immobilised coordinate frame z, a general mass balance equation for the particulate or dissolved tracer can be expressed as follows (Boudreau, 1997): @jC @ðjCUÞ @ @C þ ¼ . (A.1) jn @t @z @z @z Here C is either the volumetric concentration of the particulate tracer or concentration of the dissolved tracer in pore water, U is the material velocity given in a timeinvariant coordinate frame, and j is the porosity when C is a dissolved tracer and j1 when (1) is applied to a particulate fraction. Introducing a new set of independent variables t ¼ t; z ¼ z ðtÞ;
@z ¼1 @z
(A.2)
and applying a chain rule to link derivatives in the old system to those in the new system @jC @jC @jC @z @jC @jC ¼ ¼ ; þ @t @t @z @t @z @z the mass balance Eq. (A.1) is transformed into @jC @jC @z @jCU @ @C þ ¼ jn . þ @t @z @t @z @z @z
(A.3)
(A.5)
(A.6)
and suppressing the star superscripts in the notation of the independent variables, equation for the mass conservation in a time-varying coordinate frame can be expressed as @jC @ @ @V g ¼ C jC . (A.7) j U˜ n @t @z @z @z Here U˜ ¼ U V g is the material velocity expressed in a time-varying frame. Assuming C ¼ 1, and zero diffusion coefficient in (A.7) gives a mass conservation equation for the liquid phase w
@j @ðjU˜ Þ @V g þj þ ¼ 0. @z @t @z Equation for the concentration of solids reads @V g @C @ @ ¼ U˜ n C C . @t @z @z @z
Introducing a time-varying coordinate frame Vg ¼ U, and substituting in (A.10) porosity by void ratio (1j) ¼ (1/(1+e)), reduces (A.10) to @ @U . ¼ ð1 þ Þ @z @t
(A.11)
Using the expression for the sediment velocity in a cohesive consolidating sediment (Toorman, 1996) K @se @ K h ðrs rw Þ U ¼ hw (A.12) gr @ @z ð1 þ Þ rw and introducing new coordinates @x/@z ¼ (1/(1+e)), transforms (A.11) to the finite strain consolidation equation derived by Gibson et al. (1967) @ @ Kh @se @ K h ðrs rw Þ ¼ . w w @t @x ð1 þ Þg r @ @x ð1 þ Þ r (A.13)
(A.4)
Introducing the velocity of new coordinate levels as seen by an observer in an immobilised coordinate system (note that @z*/@t is velocity of the old frame points with regard to new coordinate levels) @z , @t
side of (A.9) by the sediment density, and neglecting diffusion, transforms (A.9) into
@V g @ð1 jÞ @ ð1 jÞU˜ þ þ ð1 jÞ ¼ 0. (A.10) @t @z @z
Here se is an effective stress, and Kh is hydraulic conductivity.
After a number of manipulations involving expansion of partial derivatives and combining different terms, (A.4) can be transformed to @jC @ @z @ @ @z þ U n C jC ¼ 0. j þ @t @z @z @z @t @t
Vg ¼
1175
(A.8)
(A.9)
Having sediment represented by single class of particles: (C/r) ¼ (1j), and dividing the left and right hand
References Ariathurai, R., Krone, R.B., 1976. Finite element model for cohesive sediment transport. Journal of the Hydraulics Division 102, 323–338. Blumberg, A., Mellor, G., 1983. Diagnostic and prognostic numerical circulation studies of the South Atlantic Bight. Journal of Geophysical Research 88, 4579–4592. Boudreau, B.P., 1997. Diagenetic Models and Their Implementation: Modelling Transport and Reactions in Aquatic Sediments. Springer, Berlin, 414pp. Dyer, K.R., 1989. Sediment processes in estuaries: future research requirements. Journal of Geophysical Research 94 (C10), 14327–14339. Gibson, R.E., England, G.L., Hussey, M.L., 1967. The theory of onedimensional consolidation of saturated clays-I. Geotechnique 17, 261–273. Govindaraju, R., Ramireddygari, R., Shrestha, P., Roig, L., 1999. Continuum bed model for estuarine sediments based on nonlinear consolidation theory. Journal of Hydraulic Engineering 125, 300–304. Harris, C., Wiberg, P., 2001. A two-dimensional, time-dependent model of suspended sediment transport and bed reworking for continental shelves. Computers & Geosciences 27, 675–690. Herzfeld, M., 2006. An alternative coordinate system for solving finite difference ocean models. Ocean Modelling 14, 174–196. Luff, R., Moll, A., 2004. Seasonal dynamics of the North Sea sediments using a three-dimensional coupled sediment–water model system. Continental Shelf Research 24, 1099–1127. Madsen, O.S., 1994. Spectral wave-current bottom boundary layer flows. In: Proceedings 24 International Conference on Coastal Engineering. American Society of Civil Engineers, Kobe, Japan, pp. 384–398. Margvelashvili, N., Saint-Cast, F., Condie, S., 2008. Numerical modelling of the suspended sediment transport in Torres Strait. Continental Shelf Research 28, 2241–2256. Mellor, G., Yamada, T., 1982. Development of a turbulence closure model for geophysical fluid problems. Reviews of Geophysics and Space Physics 20, 851–875. Mitchener, H., Torfs, H., 1996. Erosion of mud/sand mixtures. Coastal Engineering 29, 1–25. Narasimhan, M., 1992. Principles of Continuum Mechanics. Wiley, New York, 584pp. Partheniades, E., 1965. Erosion and deposition of cohesive soils. Journal of the Hydraulic Division 91, 105–139. Phillips, N.A., 1957. A coordinate system having some special advantages for numerical forecasting. Journal of Meteorology 14, 184–185.
ARTICLE IN PRESS 1176
N. Margvelashvili / Computers & Geosciences 35 (2009) 1167–1176
Reed, C.W., Neidiroda, A.W., Swift, J.P., 1999. Modelling sediment entrainment and transport processes limited by bed armoring. Marine Geology 154, 143–154. Robson, B., Webster, I., Margvelashvili, N., Herzfeld, M., 2006. Scenario modelling: simulating the downstream effects of changes in catchment land use. Cooperative Research Centre for Coastal Zone, Estuary and Water Way Management, Technical Report No.41, Indooroopilly, Australia, 34pp. Skogen, M., Svendsen, E., Bernsten, J., Aksnes, D., Ulvestad, K., 1995. Modelling the primary production in the North Sea using a coupled 3-dimensional physical-chemical-biological ocean model. Estuarine, Coastal and Shelf Sciences 41, 545–565. Smith, J.D., McLean, S.R., 1977. Spatially averaged flow over a wavy bed. Journal of Geophysical Research 82, 1735–1746.
Thomas, D.N., Judd, S.J., Fawcett, N., 1999. Flocculation modelling: a review. Water Resources 33, 1579–1592. Toorman, E., 1996. Sedimentation and self weight consolidation: general unifying theory. Geotechnique 46, 103–113. Torfs, H., Mitchener, H., Huysentruy, H., Toorman, E., 1996. Settling and consolidation of mud/sand mixtures. Coastal Engineering 29, 27–45. Wainright, S., Hopkinson, C., 1997. Effects of sediment resuspension on organic matter processing in coastal environments: a simulation model. Journal of Marine Systems 11, 353–368. Walker, S.J., Sherwood, C.R., 1997. A transport model of Port Phillip Bay. CSIRO Division of Oceanography, Technical Report No. 39, Hobart, Australia, 59pp. Van Rijn, L., 1993. Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas. Aqua Publications, Amsterdam (Multiple pagination).