A NOVEL CELLULAR AUTOMATA APPROACH TO ESTUARINE MORPHODYNAMIC MODELLING

The need for reliable modelling techniques for the pr diction of morphological change in coastal setti ngs has become increasingly important to coastal planners and poli cy makers in recent years due to the effects of acc elerated sea level rise and a shift in approach from coastal defence t o coastal management. In this research a new cellu lar a tomata based model is developed in order to bridge the gap between current bottom-up, process based models an d top-down behaviour oriented models of estuarine morphodynami c evolution, and make predictions of morphological hange over medium time-scales (one year to several decade s). The key processes of tidal flows, waves, sedim nt transport and salt marsh ecology are represented in simplifie d form in order to capture the complex interactions a d feedbacks that occur between them and which ultimately determ ine how the morphology will develop in response to environmental change. The initial bathymetry of the estuary is represented by a cluster of rectangular cells in a CA domain. Tidal flows are estimated using a new rout ing model, which shows good agreement with a conven tional 2D hydrodynamic model but is much more computationally efficient.


INTRODUCTION
Estuaries are frequently home to significant populations, as well as economic and recreational activity but are subject to morphological change caused by environmental changes such as sea level rise and infrastructure development (HR Wallingford, 1997).Esturine morphology is controlled by complex interactions between tides, waves, fluvial flows, sediment supply and underlying geology.It may also be significantly affected by the action of biological organisms (e.g.salt marsh).Robust modelling tools are needed to aid the management of risks such as flooding and threats to habitats in these areas.
Current modelling techniques include process based "bottom-up" modelling and "top-down" methods, which typically use data analysis and simple models based on equilibrium assumptions.Bottom-up models typically use two or three dimensional hydrodynamic models, coupled with sediment transport models and are best suited to making quantified predictions of change over small time and space scales, while top down models are used by geomorphologists to make qualitative predictions over much longer space and timescales (EMPHASYS Consortium, 2000).More recently, a number of hybrid models have also been developed, which combine elements of bottom-up and top-down models to bridge the gap in capability between these model types (Huthnance et. al., 2007).
In this paper a new Cellular Automata (CA) based model is presented that is capable of predicting morphological changes in estuaries over decadal timescales.

MODEL STRUCTURE
The model utilises simplified representations of tidal flows, waves, sediment transport and salt marsh ecology, which are inter-related as shown in figure 1.The model domain is divided into a regular grid of square cells.Tidal flows (depth and velocity), waves (height and period) and salt marsh (biomass density) are calculated for each cell, based on the output from the previous time step.Sediment transport rates are then calculated and used to update the bed level and composition at each cell.

TIDAL FLOWS
Cellular models of river braiding and landscape evolution (e.g.Murray &Paola, 1994 andCoulthard et al., 1999) have used routing schemes based on local slopes to estimate depth and velocity of flow at each cell.The advantage of these schemes over more sophisticated models is that they offer much greater computational efficiency, enabling simulations to be carried out over larger temporal and spatial scales.Since slope based routing is not suitable for estuarine tidal flows, an alternative approach based on bed friction has been developed for this model and was found to offer similar advantages.
Tidal water levels in the cellular model domain are first determined for a complete tidal cycle using a one dimensional hydrodynamic model.These are assumed to be insensitive to small changes in bed morphology but can be recalculated at suitable intervals or when cross section geometry changes have exceeded a predefined threshold value.Outflow from each cell is then determined from mass continuity (equation 1).

( )
Where Q o is flow out of a given cell, Q i is flow into that cell, A is the plan area of the cell, ∆L is the increase in water level and ∆T is the length of the time step.
The total outflow from each row of cells in the cellular domain is summed to give the total inflow to the next row.A target inflow to each cell in the next row of cells is calculated using equation (2), which is derived by equating the gravitational force with the bed friction, as determined by the Chezy equation.

( )
Where Q Tg is the target inflow to cell n in row (r + 1), m is number of submerged cells in row (r + 1), Q r is outflow from cells in row (r), h is water depth for cells in row (r+1) and C is the Chezy coefficient.
Transverse flows are applied within the current row as required to achieve the target distribution of axial flow in the next row, subject to constraints based on the water depth and axial flow rate.The maximum transverse flow from a cell is calculated using equation (3).

(
) Where F T is a model parameter, Q i is the axial inflow and Q i(Tg) is the target inflow.Q s is an additional allowance based on the change in storage at the current row, so that small flows over the tidal flats can occur at larger angles and is based on the change in storage within a single row of cells.The value of F T is typically between 0.25 and 1.0; larger values may give better results in areas where larger transverse flows would normally occur (e.g.where flow divergence occurs near the mouth of an estuary).However since the model does not allow for the additional frictional drag due to transverse flows, larger values can result in unwanted unrealistic morphological feedback effects (such as zigzagging).
Since the transverse flow constraints can result in axial flow inputs to the next row that are larger than the corresponding target flow calculated from equation ( 2), it is necessary to also specify limits to their magnitude.At present these have been defined in terms of velocity and the target flow in the next row (equation 4).
Where F A is a constant model parameter, Q Tg is the target flow, Q i is the calculated inflow in the current row, H (r) and H (r+1) are the water depths in the current and next rows and w is the cell width.The value of F A is between 1 and about 1.5.A larger value gives the model more freedom to adjust to changes in bed profile without carrying out the additional computation required when Q A(Max) is exceeded.
If the limits on axial flow magnitude are exceeded then the excess flow is redistributed into the cells with additional axial flow capacity.The redistributed flows are routed back upstream along with the negative residual flow in the cells which had exceeded the axial flow limits.This process allows the pattern of flow to be adjusted upstream of any sudden change in bathymetry, as shown in figure 2.
Model performance was tested by comparing results with those from TUFLOW (Syme, 1991), a commercial 2D hydrodynamic model.Some example comparisons are presented in figure 3.

WAVES
Wave height and period in the basin of the estuary are estimated using the wave hincasting equations given in the Coastal Engineering Manual (US Army Corps of Engineers, 2002).Fetch is estimated at each cell using the procedure described below, while wind speed and direction are supplied as an input at each time step or generated using a simple statistical wind model.Only locally generated waves are considered by the model.Fetch is estimated using an accounting procedure, starting from the 'up wind' model boundary.This is the first or last row or column depending on the wind direction as indicated in figure 4. For wind directions between -45 and +45 degrees the procedure begins with row n, where each cell is assigned a boundary fetch (a fixed model parameter).This fetch is reduced if necessary to limit the resulting wave height to the depth limited wave height (taken to be 0.78H).By reducing the fetch, as opposed to simply depth-limiting the calculated wave height, the wave heights for any 'down wind' cells with larger depths are also reduced.The fetch from each cell is then distributed to seven cells in row (n -1), in proportion to the subtended angles; this process smoothes the results by averaging and allows an approximate representation of diffraction around sand bars and other obstacles.An offset is applied periodically in this process to account for the difference between the grid and wave directions (e.g. for a wave direction of 45 degrees the offset is applied at every row and for a direction of 26.5 degrees it is applied approximately every second row).The wave direction is initially assumed to equal to the wind direction.
The fetch for cells in row (n -1) is determined from the sum of the contributions from cells in row n plus the additional fetch (the length of one cell (L) for a wind direction of zero degrees).The wave direction is determined from the average of the directions of the contributions of fetch to that cell, weighted according to the magnitude of each contribution.The procedure is repeated to obtain the fetch for cells in row (n -2) and subsequent rows.Wave height and period are then calculated for each cell using equations 5a and 5b (US Army Corps of Engineers, 2002).
Where H m0 is the energy based significant wave height, T p is the wave period, X is the fetch and u * is the friction velocity

SEDIMENT TRANSPORT
The Van Rijn TR2004 formula (Van Rijn, 2007a, b, c) is used to calculate sediment flux at the centre of each cell in the model domain.This method is capable of determining transport rates of single or multi-fraction sediment due to both waves and currents.Bed load is calculated using equation 6a.
Where γ is a coefficient (0.5), ρ s is the sediment density, d 50 is the median particle size, D * is a dimensionless particle size, τ' b,cw is the instantaneous bed shear stress due to both waves and currents and τ b,cr is the critical bed shear stress.Net bed load transport is obtained by time averaging over the wave period.For suspended load calculation the reference concentration is given by equation 6b.
( ) Where p clay is the proportion of clay in the sediment and f silt = d sand /d 50 , d sand = 62µm.T is the dimensionless bed shear stress parameter (which is related to the bed shear stress and critical bed shear stress) and a is the reference level (which is related to the bed roughness).Suspended sediment transport is then calculated by integrating the concentration and velocity profiles over the depth.
When operating in single fraction mode, bed levels are updated in accordance with equation 7 (Soulsby, 1997).
Where z is the bed level, q tx and q ty are the total volumetric sediment transport rates and ε is the bed porosity.Sediment is moved into neighbouring cells in proportion to resolved components of the flow vectors, as shown in figure 5. Additional sediment transport along lateral slopes is included and is assumed to be proportional to the difference in bed level and the principle transport rate.This is similar to the lateral transport rule used in previous river braiding cellular automata models (e.g.Murray and Paola, 1994).
For multi-fractional sediment transport modelling, sediment is moved into neighbouring cells in accordance with the calculated bed load and suspended load transport rates for each fraction.Erosion and deposition are calculated separately for each fraction based on the difference between the amount of sediment transported into and out of each cell.Bed composition is stored in a series of sub-layers and a top 'active' layer.Deposited or eroded sediment is added to or removed from the active layer as appropriate, however in the case of deposited suspended sediment it is first added to a deposition layer from which it is transferred to the active layer at a rate determined to the fall velocity and diffused into the neighbouring cells (advection is not represented for the deposition layer because this would require much shorter time steps to satisfy the Courant number criterion).The active layer is restored to its default thickness at the end of each time step by adding or removing material to or from the sub-layers.

SALT MARSH DYNAMICS
Salt marsh has a significant impact on the morphological evolution of many estuaries.Marshes stabilise and trap sediment but are susceptible to erosion and drowning due to changes in sea level and changes to wave climate.A simple representation has been included in the model in an effort to capture some of the interactions and feedbacks that commonly occur between estuarine morphology and salt marshes.
Marsh biomass is estimated from the proportion of time for which a cell is inundated during the previous month.Mariotti and Fagherazzi (2010) related biomass to the difference between marsh elevation and highest astronomical tide; however since the model must allow for changing tidal levels, due to sea level rise and/or tidal propagation effects, the equivalent measure of inundation time (T i ) was chosen (Mudd et al., 2004).Mariotti and Fagherazzi (2010) proposed a parabolic relationship, with biomass falling to zero at both the maximum and minimum inundation levels; however Mudd et al. (2004) propose a linear function with biomass increasing from zero at the minimum inundation level (T i(min) ) to a maximum at the maximum inundation level (T i(max) ), with biomass then falling to zero for inundation greater than the maximum.D' Alpaos et al. ( 2006) have suggested a similar linear function for cases where the marsh is dominated by a single vegetation species (Spartina alterniflora); however they find that in cases where a variety of halophytic species are present, a different relationship exists.In these cases they propose a linear relationship with biomass increasing from zero at T i(max) to a maximum value at T i(min) and above.This latter relationship has been initially adopted for the CA model.Mudd et al. (2004) have proposed a maximum value of 0.6 for T i(max) , however in practice this value could be expected to vary according to local conditions and it has therefore been left as an adjustable parameter in the current CA model setup.
The biomass density is used to calculate the increase in critical bed shear stress, enhanced sediment trapping and organogenic biomass generation using the relationships given by Mariotti and Fagherazzi (2010).The relationship between biomass and hydraulic roughness is assumed to be linear, with Chezy C values between 50 (zero biomass) and 10 (maximum biomass).These values were estimated with reference to Manning's n values (using the approximation C = 1/n) for similar conditions, as given by Chow (1959).

SENSITIVITY TESTING
A number of sensitivity tests have been carried out to assess the affects of the various model parameters and components.A simplified estuary bathymetry was created for this purpose, with dimensions 15km by 1km and divided into 6000, 50 by 50m, grid cells, arranged in 300 rows and 20 columns.A tidal flat set at approximately neap high tide level was incorporated together with a single channel of varying width and depth.The Model was run to simulate a period of 10 years of morphological evolution with fluvial inflow, tidal variations, and sediment properties varied individually in separate tests.The tests were initially carried out using a single fraction sediment transport approach but were repeated using a multi-fraction approach and with and without the salt marsh model.Cross sections showing preliminary sensitivity test results are shown in Figures 7a to 7d.These figures show the differences in the cross section profile after 10 years simulation time, when fluvial X inflow or tidal range is increased.In the base case scenario, the fluvial inflow is 1m 3 /s and the tidal range is between 2m (neap) and 4m (spring).The fluvial inflow is then increased to 10m 3 /s and in a separate test the tidal range is increased to between 3m and 6m.In all scenarios the mean sea level is set at 49m (arbitrary datum).
Results for the single fraction and multi-fraction model versions are given separately.The median sediment size (D 50 ) is 0.2mm in all cases and D 90 is set at 0.5mm.For the multi-fraction model sediment is divided into three fractions: Fraction 1: 0.10 -0.15mm (33.3%)Fraction 2: 0.15 -0.25mm (33.4%)Fraction 3: 0.25 -0.55mm (33.3%) Results for the multi-fraction model simulation are broadly similar to those for the single fraction version.In both versions small amounts of deposition has occurred on the tidal flats and this was slightly higher for the multi-fraction model version, which is presumably due to the higher mobility of the fine sediment fraction.It was noted that the median sediment size had increased in areas of erosion and reduced in areas of deposition as might be expected.This is shown in Figure 8.
Although the difference between the results for the multi-fraction model is small for the current set of tests, it is expected that where silt and very fine sand particle sizes are present the multi-fraction capability will be more significant and will enable the simulation of sediment sorting and deposition of fine particles onto the tidal flats.transport components have been successfully implemented.The model is computationally efficient, which will be a very useful feature in long term simulations.
The next stage of the model development will be to investigate rules for sediment flux at the downstream boundary so that simulations can be made to evolve towards a state of dynamic equilibrium.It should then be possible to test the response to changes, such as the addition or removal of sediment or an increase in mean sea level.

Figure
Figure 1: Model Structure

Figure 2 :
Figure 2: Tidal Model Example Results

Figure 5 :
Figure 5: Sediment transport directions within the CA

Figure 6 :
Figure 6: Contour map showing initial conditions (left) and base case results (single fraction method).The approximate location of the time series flow results in Figure 3 is indicated by X.The dashed lines show approximate positions of the cross sections given in Figure 7.

Figure 7d :Figure 8 :
Figure 7a: Single Fraction Sensitivity Test Rests (Inflow), at chainage 500m from upstream extent of model domain