Modelling the sensitivity of suspended sediment profiles to tidal current and wave conditions

Seawater turbidity due to suspended particulate material (SPM) is an important property of a marine ecosystem, determining the underwater light environment and many aspects of biological production and ecology. SPM concentrations are largely determined by patterns of sediment resuspension from the seabed due to shear stress caused by waves and currents. Hence planning for the construction of large scale offshore structures which will alter regional hydrodynamics needs to consider the consequences for SPM concentrations. Here we develop a one-dimensional (vertical) model of SPM dynamics which can be used to scope the effects of changes in wave and tidal current properties at a site. We implement the model for a number of sites off the east coast of Scotland where we have extensive data sets to enable numerical parameter optimisation. The model performs well at simulating fluctuations in turbidity varying from flood-ebb tidal cycles, spring-neap cycles, storm wave events, and an annual cycle of SPM concentration which is attributed to seasonal consolidation of seabed sediments. Sensitivity analysis shows that, for the range of seabed sediment types in the study (water depth 16e50 m; mud content 0.006e0.380 proportion by weight), relatively large (50%) attenuations of tidal current speed are required to produce changes in water column turbidity which would be detectable by observations given the variability in measurements. The model has potential for application to map the large scale sensitivity of turbidity distributions to the installation of wave and tidal energy extraction arrays. © 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Sea water turbidity due to suspended particulate matter (SPM) determines the depth to which sunlight penetrates below the sea surface. This is one of the key factors determining the species composition and productivity of marine ecosystems. The effects include the rate and fate of primary production, the performance of visual predators such as fish, potential for refuge from predators by vertically migrating species, and the scope for seabed stabilisation by algal mats. Hence, turbidity is a key property of an ecosystem, but one which has proved to be particularly difficult to model in shelf and coastal systems.
Some of the material contributing to turbidity may be of biological origin, but in coastal waters the majority is mineral particles originating ultimately from seabed disturbance and land erosion, the latter being deposited in the sea by rivers and aerial processes. SPM is maintained in the water column or deposited on the seabed depending on combinations of hydrodynamic processes including baroclinic (density-driven) or barotropic (mainly tidal and wind driven) currents, and wave action (Ward et al., 1984;Huettel et al., 1996). Spatial and temporal variations in hydrodynamics, or interventions such as engineering structures which alter hydrodynamics, should therefore be a major determinant of turbidity.
Full simulation of the impact of waves and currents on suspended sediment concentrations requires the solution of equations representing erosion and deposition of sediment from the seabed, together with vertical mixing and horizontal transport in the water column. Typically the mixing and advection terms are posed as partial differential equations embedded in a computational scheme for solving the equations of fluid dynamics (e.g. Teisson, 1991). There are several systems available for this task (e.g. Gerritsen et al.,, 2000;Mercier and Delhez, 2007;Warner et al. 2008;Danish Hydraulics Institute, 2013). However, in each case the inclusion of SPM simulation adds considerably to the computational demand and requires extensive and costly calibration of area-specific parameters. For many applications, this may be prohibitively demanding. Some authors have explored alternative 'short-cut' approaches involving e.g. blending of satellite remote sensing data on SMP concentrations and simulated hydrodynamic flow fields (Wu et al., 2011). Here, we propose a 'lightweight', onedimensional (vertical), modelling approach for basic simulation of SPM dynamics, incorporating simple caricatures of the fundamental erosion and deposition processes which can be used to quickly scope the effects of hydrodynamics on turbidity distributions. Our approach is to simulate time-dependent vertical profiles of suspended sediment concentrations at point locations, given seabed depth and mud content, and time-dependent bed shear stress and sediment erodibility. Clearly, this approach cannot take account of lateral transport of suspended sediment, so its use must be limited to areas where the majority of sediment material in the water column arises from local seabed resuspension rather than horizontal transport.

Key processes affecting the vertical distribution of suspended sediment
In a closed, one-dimensional (vertical) system the mass of SPM in the water column represents the balance between erosion and suspension rates of seabed sediment, and deposition rates of suspended material. The main proximate drivers of these rates are timevarying vertical diffusivity and shear stress arising from friction between the seabed and flowing water, in particular the orbital flows which occur beneath surface waves, and directed flows due to tides and residual currents. However, the context is set by a variety of seabed sediment properties including bedform architecture, grain size composition, cohesion, consolidation and compaction. Cohesion arises primarily from electrochemical attraction forces between particles, compaction from gravitational compression leading to extrusion of pore-waters, and consolidation from adhesion forces between particles due to inorganic chemical reactions and organic molecules produced by microbiological activity. In addition, bioturbation of sediments by sifting and burrowing fauna may lead to modification of erodibility.
The shear stress on a seabed particle is a function of its size, the flow speed, and the densities of the fluid and particles (Wilcock et al., 2009). When the shear force exceeds resisting forces due to gravity, cohesion and consolidation, then a particle can become mobile. As shear forces increase, particles initially undertake short hops along the seabed (saltations), or rolling motions. Such particles are said to be part of the 'bed-load'. When the value of the bed-shear velocity becomes sufficiently high relative to the particle fall velocity, then bedload particles can be lifted into suspension. The vertical flux of particulate mass can be described by the differential equation: where C(z) is the suspended sediment concentration at altitude z above the seabed, C a is the concentration at a reference altitude z a close to the seabed, K s is the vertical diffusivity, and ω s is the fall velocity of particles. Predictions of vertical distributions of concentration therefore depend on assumptions about the vertical profile of diffusivity. Commonly used alternatives are to assume a constant diffusivity with altitude above the seabed, a linear increase, or a parabolic variation with peak diffusivity in mid-water. Assuming a linear increase with altitude, the concentration profile is given by where u* is the shear velocity at the seabed, κ is the von Kármán constant (0.4), and β is a coefficient relating eddy viscosity to eddy diffusivity (taken to be 1) (Rouse, 1937;Van Rijn, 1984;1993). The exponent ω s /(β·κ·u*) is referred to as the Rouse number. Alternative assumptions regarding the vertical distribution of diffusivity give different expectations for the vertical profile of concentration, but the linear Rouse approach is most commonly applied (Camenen and Larson, 2007).
Sinking velocity is a critical term for both the initiation of particle motion on the seabed, and the structure of vertical profiles of SPM concentration in the water column. At equilibriumwhere the sum of the gravity force, buoyancy force and fluid drag force are equal to zero -the downward sinking velocity of particles depends on the density and viscosity of the fluid, and the density, size, shape, and surface texture of the particle. The classical Stokes equation for the fall velocity of a particle assumes a spherical shape and laminar flow (Reynolds numbers less than 1). Despite extensive research there is still no analytical solution to predict the fall velocity of natural shaped particle, or particles large enough to generate turbulent flow (Camenen and Larson, 2007). Many investigators have proposed empirically based relationships to predict particle fall velocities with varying degrees of complication and success (Sadat-Helbar et al., 2009).
Although particle shape is certainly a factor contributing to uncertainty in sinking rates, part of the variability arises from particle-particle collisions during suspension in the water column. Collisions of fine grained particles can lead to aggregation and formation of flocs with potentially enhanced sinking rates, depending on the physical cohesive properties of particle grains and their stickiness due to biological coatings (e.g. Krone, 1978;Mehta, 1989;Andersen and Pejrup, 2002;Winterwerp, 2002;You, 2004). The probability of collisions will be a function of the suspended sediment concentration. Experimental studies have found that settling velocity for mud and silt particles is independent of concentration below 0.4 g/l. Between 0.4 and 2.0 g/l, settling velocity increases with concentration due to flocculation. Above 2.0 g/l settling velocity rapidly decreases due to the break-up of flocs, mutual hindrance, and interactions between the flows around adjacent flocs that tend to increase upward friction (Cancino and Neves, 1999). A widely used empirical relationship describing this process (Burt, 1986) is of the form: where k and γ are constants, and C lies between a lower threshold for particle-particle interactions, and an upper threshold at which particles begin to interfere and the effective settling velocity is reduced. The upper concentration corresponds to values found in e.g. mud slides, where the water-sediment mixture forms a super-dense liquid which dampens turbulence and reduces shear stress as a feedback process (e.g. Richardson and Zaki, 1954). Whilst this phenomenon may occur in highly turbid estuaries, it is not expected in typical shelf-sea marine situations and we do not take it into account here.
The velocity of suspended particles in a longitudinal direction is almost equal to the fluid velocity. Lateral transport of suspended sediment is therefore simply the product of the vertical profile of sediment concentration and the vertical profile of water velocity (Van Rijn, 1993). Hence, horizontal bed-load transport in fluctuating flow regimes is relatively easily modelled because vertical processes affecting the particles are limited to the onset and cessation of motion on the seabed. However, modelling suspended loads in fluctuating flows is more complicated because deposition and erosion fluxes are decoupled in time depending on the height of the water column and the rate of vertical diffusivity. Full modelling of suspended sediment transport therefore requires dynamic representation of vertical convection-diffusion processes in order to resolve short term fluctuations in vertical concentration gradients. However, by far the most difficult aspect of the problem is representation of the seabed sediment context. The widely used Shields relationship which sets a critical shear for particle motion depends only on grain size (Shields, 1936;Paphitis, 2001;Beheshti and Ataie-Ashtiani, 2008). No general relationships have emerged to represent the effects of consolidation effect on sediment erodibility (McCave, 1984). Early formulations which simply propose a site and time specific parameter to represent erodibility (Partheniades, 1965) remain widely used in models of sediment processes (e.g. Whitehouse et al., 2000;Ribbe and Holloway, 2001;Kuhrts et al., 2004;Pandoe and Edge, 2004;van den Eynde, 2004).

One-dimensional suspended sediment model specification
The model caricatures three main processes leading to the generation of vertical profiles of turbidity:  The availability of fine-grained material on the seabed which may contribute to turbidity in the water column  Effects of fluctuations in bed shear stress on the balance between net erosion and net deposition  Balance between vertical mixing and sinking of particles suspended in the water column Parameters and terms are listed in Table 1. The model does not take account of horizontal transport processes.

Availability of fine-grained material
We assume that water column turbidity is mainly due to suspended mud and silt grains, i.e. particles <0.063 mm diameter, and that the availability of these particles for lifting off the seabed is expressed by the product of some power of the weight-specific mud content of the sediment (S ε ), and an erodibility term (E β ).
The erodibility of the mud fraction of sediment depends on a variety of factors but we assume it to depend mainly on consolidation due to biological activity. We expect this to follow a seasonal cycle dictated by temperature and the input of fresh organic matter settling from the spring and summer plankton blooms. We do not know the exact form of this, though observational data on phyto-detritus pigments in North Sea sediments, oxygen consumption and nutrient fluxes indicate a peak of activity in June/July and a minimum in December/January. In addition, we know that pigment concentrations and microbial fluxes increase with the mud content of sediments (Serpetti, 2012;Serpetti et al., 2012). So, we caricature the erodibility of sediments by a time dependent cosine function scaled to vary between arbitrary non-zero, positive limits (0.5 and 1.0), and phase shifted by a period t β relative to the annual solar cycle: where t is the time in days. The availability of fine-grained material for suspension into the water column is then given by (S ε ·E (t) β )

Fluctuations in erosion and deposition rates
Variations in bed shear velocity lead to fluctuations in the rates of erosion from the seabed and deposition from the water column. However, the key point is that vertical deposition fluxes reflect time-lagged signals of past erosion events. Explicit simulation of these processes requires the computational solution of fluid dynamics equations including the advection and diffusion of suspended particles. Here, we caricature the net effect of these processes by posing that the bed stress forming any given instantaneous vertical profile of suspended sediment is a time-weighted average of the stress over some period prior to observation.
We define an exponentially declining time-weighting function where t is a vector of shear observation times prior to the instant at which an observation is made, T a ≥ t ≥ 0, and T a is a negative number representing the autocorrelation time scale relevant to the formation of the suspended sediment profile.
The time-weighted shear (τ a(t) ) is then given by where τ (t) is a vector of bed shear stress at time t.
The corresponding time weighted bed shear velocity is then given by: where ρ is the density of seawater.
Then, we represent the near-bed suspended sediment concentration by: This expression contains three components: the scaling coefficient α which equates the modelled concentration to observed measurement units; the term for availability of finegrained material (S ε ·E (t) β ), and the time-weighted bed shear stress term (τ a(t) δ ) which corresponds to the erosion rate expression of e.g. Partheniades (1965). We do not set an explicit threshold of shear stress for the initiation of particle motion, since we are not addressing sediment erosion fluxes or steady state concentrations under constant flows. Rather, we aim to caricature transient concentrations in a time varying system, where the concentration near the seabed at any instant reflects the balance between deposition and erosion fluxes.

Balance between vertical mixing and sinking of suspended particles
We assume a linear profile of vertical diffusivity, and hence that the distribution of suspended material can be primarily explained by the Rouse formulation. Hence, the suspended sediment concentration at altitude z (0 < z ≤ h), is given by: The exponent here corresponds to the Rouse number, but additionally incorporating the expression ( ( ) ) to reflect increasing particle aggregation in the water column with increasing sediment concentration (Burt, 1986).

Model parameter fitting and validation
The model contained 9 parameters that required to be estimated by fitting to observed data ( Table 1). As an observed dataset, we assembled measurements of turbidity at 0.5 m depth intervals from 10 sites off the east coast of Scotland, sampled at varying intervals between mid-2008 and the end of 2011. The data set was divided into two approximately equal subsets in terms of number of observations, on the basis of sampling date. The earlier subset was treated as calibration data to which the model was fitted. The later subset was treated as validation data.
All 9 parameters of the model were fitted by minimising the r.m.s error between the entire calibration set of observed turbidity-at-depth, and predicted values assuming the inputs of bed shear stress time series, seabed mud content, and sea surface altitude above the seabed at each site. Minimisation was performed by standard Nelder Mead optimisation using the 'optim' function in the R statistical environment, with hessian matrix output so as to derive the standard errors of the parameters. The quality of the fit was measured with the Pearson correlation coefficient.
The fitted parameters of the model were then used to predict the time series of turbidity at two horizons in the water column at each site (5 m altitude above the seabed, and 5 m depth below the sea surface) for the full duration of the sampling period. The predictions for the calibration and validation period at each site where then compared with the measured turbidity using the Pearson correlation coefficient.

Validation of fitted sinking rate
Particle sinking rate was the only fitted parameter of the model which could be independently validated from empirical evidence. We did so by estimating the particle diameter of the suspended material in the water column implied by the fitted sinking rate, assuming two alternative empirically-based relationships between sinking rate and diameter. We expected the implied particle diameter to be less than the 0.063mm threshold for mud grains.
Sadat-Helbar et al. (2009) reviewed 17 published relationships between sinking rate and particle diameter and identified the formulation developed by Wu and Wang (2006) as being one of the most realistic: where A, B and a are coefficients, d is the actual particle diameter, the term D gr is referred to at the effective grain size, and v is the kinematic viscosity. Empirical calibration against a wide range of sediments provided coefficient values as: where S f is the Corey shape factortypically taken to be 0.7 (Camenen, 2007), and = • ( Here, is the density of the sediment material, is the fluid density, and g is the acceleration due to gravity.. -Helbar et al. (2009) also provided their own generalised piecewise relationship in which fall velocity increases as a power function of particle diameter, without incorporating any shape parameter terms:

Study region and data for application of the model
Data to drive, calibrate and validate the model came from inshore waters (up to 12 km from the coast) off Stonehaven on the North Sea coast of Scotland ( Figure 1). Turbidity and seabed sediment data came from field observations at a set of 10 sites which were sampled at various frequencies over the period July 2008 -December 2011. Data on wave and current properties came from a hydrodynamic model of the region.

Background data on the sampling sites
Data on the bathymetry and sediments of the study region have been published elsewhere (Serpetti, 2012;Serpetti et al., 2012). Very briefly, detailed bathymetric data were collected during December 2006 using a Simrad EM950 95 kHz multibeam sonar (Serpetti et al., 2011). Seabed sediment sampling for grain size and chemical analysis was subsequently carried out at >50 sites in this region using a 0.1 m 2 Day grab during two surveys in April 2007 (RV 'Clupea') and September 2008 (RV 'Alba na Mara') 2012). Of these, 7 sites were selected for approximately monthly measurements of water column turbidity and sediment porosity, permeability, oxygen consumption and nutrient fluxes between mid-2008 and mid-2009 from the RV 'Temora'. At these sites, undisturbed sediment cores were collected with a Mini Muc k/MT 410 corer fitted with 60 cm acrylic core tubes (Serpetti, 2012). Water column turbidity data collected at these 7 sites, plus 2 others ("main" and "inner") which were visited weekly throughout 2008-2011, and an anchor station at which turbidity was monitored at high temporal resolution over 4 overnight periods in 2008, were used for our study reported here (Table 2, Figure 2). Sediment mud content (proportion by weight of grains <0.063 mm) at the 10 sampling sites varied between 0.006 to 0.380 and the water depth between 16 and 50m (Table 2).

Measurement of turbidity data
On each sampling occasion, turbidity (Formazine Turbidity Units (FTU), proportional to SPM (g.m -3 )) was measured at 0.5 m depth intervals between the sea surface and around 1m above the seabed using a Saiv SD 204 CTD unit (Saiv A/S Environmental Sensors & Systems) fitted with an optical backscatter sensor. Weekly sampling with this equipment was carried out from a small research vessel operated by Marine Scotland Science at the "main" and "inner" sites which were up to 5 km directly offshore from Stonehaven (Table 2, Figure  2,). One of these sites ('main' site) has been monitored weekly since 1997 (Bresnan et al. 2008). The other ('inner' site) has been sampled only since 2007. The same methodology was used to collect turbidity profile data during an intensive study between 25th and 30th September 2008 by RV 'Alba na Mara', and during each approximately monthly visit to the 7 sediment core sampling sites in 2008 and 2009 by RV 'Temora' (Serpetti, 2012).
During the intensive sampling in September 2008 from RV 'Alba na Mara', the vessel anchored each night in Stonehaven Bay. The anchoring position was not precisely the same each night, but was within an area of 400 x 400 m, with water depth 14-19 m. Throughout the periods at anchor data from a range of sensors, including a 25 cm path-length Seatech transmissometer, fed with water from the vessel's pumped seawater supply (intake depth 3m), were recorded at 1 minute intervals. The transmissometer data were subsequently calibrated in terms of beam attenuation and then re-scaled to turbidity units by intercalibration with data from the Saiv CTD system. The time series of turbidity at the anchor stations were then used as a further validation of the model.
Temperature and salinity data recorded by the system were calibrated from reversing oceanographic thermometer readings and salinometer analyses of water samples collected from near-surface and near-seabed depths on each visit to the 'main' sampling site.

Bed shear stress (τ t ) estimated from time series of modelled and observed tidal current and wave properties
Time series of bed shear stress required as input to the model could be acquired by various means: for example, computed from direct observations of current velocity and wave properties, computed from modelled tidal harmonics and wave-field hindcasts, or extracted from temporally-explicit hydrodynamic model simulations. Here, we used a combination of modelled tidal harmonics, direct observations and wave hindcasts, to generate time series of depth averaged tidal current and seabed wave orbital velocities at each of the turbidity sampling sites. Bed shear stress was then computed from the combination of tidal and wave orbital velocities. A disadvantage of this approach is that we could not account for shear arising from wind or buoyancy-driven residual currents. On the other hand, the main advantage was that we were not constrained by having access to full hydrodynamic and wave model simulation outputs for the entire duration of our study period. Using our method we had the capability to construct tide and wave-driven shear stress series for any period during which wave-monitoring buoy data were available.
We used a calibrated, 3-dimensional, coupled wave-current model for the region constructed in MIKE by DHI (MIKE 3 FM for tidal and wind-driven currents and MIKE 21 SW for spectral waves; Sabatino et al., 2016) to simulate current speed and direction, and wave parameters during a 7-month period in 2010. The model was based on an un-structured grid which varied in spatial resolution from <100 m in the interior of the study region, to >5 km on the far-field regions more than 100 km distant ( Figure 2). The MIKE 3 FM model was forced at the boundaries by sea surface elevations from the Oregon State University Tidal Prediction Software (OTPS; Egbert et al., 2010), and meterological data from the ERA-Interim analysis of wind and sea-level pressure (Dee et al., 2011). The MIKE 21 SW model was forced by wind velocity data as for MIKE 3, and swell wave conditions at the open boundaries taken from a North Atlantic scale wave model (Venugopal and Nemalidinne, 2014;. Coupling between the models was uni-directional, so that current fields affected the simulated wave spectra, but not vice-versa. The models were calibrated against seasurface elevation data from a tide gauge at Aberdeen (15 km north of the sampling sites), harmonic components of current velocities from archived recordings at 35 recording current meter (RCM) mooring locations in the region, and wave data from the UK Wavenet Firth of Forth monitoring buoy approximately 50 km from the study area (56° 11.31'N, 002° 30.43'W, www.cefas.co.uk/publications-data/wavenet/; Figure 1).
For each location of interest, we extracted the parameters of tidal harmonics from equivalent locations in the model output. These were then used to reconstruct 15 minute interval time series of depth averaged current speed for the entire period June 2008 to December 2011. We also used the model output to establish predictive linear regressions for significant wave height, peak wave period and peak wave direction at each sampling site based on the corresponding observations at the UK Wavenet Firth of Forth monitoring buoy. We then used these to reconstruct extended time series for each site using Wavenet data for the entire 2008-2011 period.
Time series of wave orbital velocities at the seabed were derived from the estimated 15 minute significant wave height and peak wave period at each site using the algorithm of Soulsby (2006). Shear velocity due to the tidal flow was calculated from the vertically averaged current speed throughout the water column using the "law-of-the-wall" method (Soulsby and Clarke, 2005) which assumes a logarithmic decrease in velocity with proximity to the sediment-water interface. Combination of tidal current and wave-induced bed shear stress was then performed according to the algorithm detailed in Soulsby and Clarke (2005) taking account of the relative directions of the currents and waves at each 15 minute interval.

Sensitivity to reductions in bed shear stress
In order to scope the impact on turbidity of reductions in tidal current speed or wave height due to energy extraction, we re-ran the bed shear stress calculation using the MIKE by DHI simulation outputs for the sampling sites but assuming some removal of either tidal power by diminishing the depth mean current speed, or wave power by diminishing the significant wave height (but not the wave period).
Provided that the water depth is larger than half the wavelength, the power associated with a wave train is where P w is the power per metre of wave front (W.m -1 ), H is the wave height and T p is the wave period. Hence, the change in wave height associated with extraction of energy so as to reduce the natural wave power by a fraction k w = P w(exploited) / P w(natural) without affecting the wave period is simply √ The equivalent measure for a current flow (power per metre at the sea surface perpendicular to the flow) is given by: where h is the seabed depth and V is the depth mean current speed. Hence the change in current speed associated with extraction of energy so as to reduce the current power by a fraction k c , is √ 3 .

Estimating the impact of changes in turbidity on light penetration depth
Prior to the study period reported here (February 2007-May 2008, vertical profiles of photosynthetically active radiation (PAR) had been collected simultaneously at the sea surface and in vertical depth profiles, during weekly visit to the main sampling site. From these data, an empirical relationship between the Beer's Law vertical attenuation coefficient of down-welling sea surface irradiation, and turbidity was established (Heath et al., 2015). The relationship also involved the in-situ concentration of phytoplankton chlorophyll (measured by a calibrated in-situ fluorometer) which absorbs a portion of the down-welling light. The fitted relationship was: where ki(z) is the light attenuation coefficient at altitude z (natural logarithmic, m -1 ), C(z) is the turbidity (FTU), and chl the phytoplankton chlorophyll concentration (mg.m -3 ) at altitude z.
Using the turbidity at 5m depth predicted by our suspended sediment model as a measure of near-surface conditions, we can use this relationship to estimate the depth at which downwelling light is attenuated to 1% of that at the sea surface, in the absence of any chlorophyll. : 1% ℎ = (0.01) 0.1473 + 0.0620 · ( ≡ 5 ℎ) 1% sea surface irradiance approximately corresponds to zero net photosynthesis i.e. gross photosynthetic uptake of carbon equals respiration. So the depth at which this occurs is a measure of the euphotic zone thickness.

Variation in wave and current properties between sampling sites
Phase and amplitude of the tidal harmonic components of depth averaged current speed and direction at locations in the MIKE model grid corresponding to the 10 turbidity sampling sites are given in Table 3. Significant wave height, mean wave period and mean wave direction simulated by the MIKE model at each sampling site, were linearly related to the temporally corresponding data at the Firth of Forth Wavenet buoy. Regression parameters for these relationships (Table 4) were used to predict the wave environment at each sampling site out-with the MIKE simulation period from the Wavenet buoy data.
Depth averaged tidal current speed was more variable between sites than the significant wave height, though the between-site differences were small in both cases. However, when combined with seabed depth and the directional data on currents and waves, the resulting bed shear stress was markedly different between sites (Table 5, Figure 3).. There was a clear relationship between the median and 95 th centile of shear stress, and the mud content of the seabed sediments (Figure 4). This leads us to conclude that the MIKE model of tides and waves, and the subsequent derivation of bed shear stress provided a realistic measure of the time variations in forces acting at the seabed.

Modelled turbidity profiles and time-series
The 371 vertical profiles of turbidity collected during the study period (30,433 individual measurements of turbidity-at-depth) were divided into two parts: data collected prior to 1 August 2009 (145 profiles, 12,044 measurements from all 9 sites, referred to as the calibration period), and data collected after 1 August 2009 (226 profiles, 18,389 measurements from the 'main' and 'inner' sites only, referred to as the validation period).
The optimised parameter set provided a statistically significant fit of the model to both the calibration and the validation data subsets. The fitted parameters and their standard errors are shown in Table 6. Figure 5 shows the scatter plot of fitted and measured turbidities (all sites, all depths) for the calibration period, and for the validation period.
The fitted sinking rate of the suspended particles (0.000210 m.s -1 , s.e. 0.000088) implied a particle diameter of 0.02 -0.03 mm which, as we anticipated, was well below the 0.063 mm upper limit for mud grains ( Figure 6). Hence we conclude that the fitted sinking rate was at least a credible value based on empirical evidence. The fitted rate equated to a period of around 3 days for a particle to sink from the sea surface to the seabed in still water conditions at most of the sampling sites, justifying the long autocorrelation time-scale of 4.7 days to emerge from the parameter optimisation.
Median and ranges (5 th and 95 th centiles) of the modelled vertical profiles of turbidity on sampling occasions during the calibration and validation periods at each site generally agreed well with the corresponding observed data (Figure 7). The model performed well at the main and inner sampling sites. High extremes (95 th centile) of observed turbidity were notably underestimated by the model at the shallowest sampling sites (core93 and core113), whilst two sites (core123 and coreRay) were insufficiently sampled to obtain a meaningful assessment of the performance of the model relative to the observed data. Figures 8 and 9 show the time series of depth-averaged tidal current speed and significant wave height that formed part of the input data to the model for the 'main' sampling site, and the fitted and observed turbidity data at two depths. Periods of extreme wave activity were under-sampled for safety reasons, but the model clearly reproduced fluctuations in observed turbidity which were associated with the spring-neap tidal cycle.

Detailed sampling at the main site during September 2008
The daily sampling at the main site during 22-30 September 2008 was carried out during a period of increasing tidal range from neap to spring tides and significant wave heights mainly less than 1 m. Earlier in the month wave heights had peaked at >4 m. Few observations were available to coincide with the period of wave activity between 4 and 14 September, but the model replicated the observed rising turbidity during 22-30 September (Figure 10).

Time-series data from the overnight anchor stations in September 2008
As at the 'main' sampling site, both the model and the observations showed an increasing trend in turbidity at the overnight anchor station in September 2008 ( Figure 11). However, the high temporal resolution observational data at the anchor site resolved 6-hourly fluctuations in turbidity associated with the flood and ebb tide which could not be resolved with daily sampling data. The model accurately replicated the timing of these fluctuations in turbidity, though the amplitude of the modelled fluctuations was generally smaller than in the observed data.

Sensitivity to seabed mud content
The sensitivity of the modelled turbidity to seabed sediment mud content is given by the term S ε . All other conditions being equal, this term implies that if the mud content of sediments at the main sampling site was 1.0 instead of the observed 0.061, then the water column turbidity would be 1.49-times higher ( Figure 12). It is not possible to visualise this response in the observed data because the wave and tidal conditions at every site are unique so we cannot isolate the effects due solely to mud content. However, the contrast in mud content between the sites implies a range of variation in turbidity of around 0.75 to 1.3-times the turbidity at the main site.

Sensitivity to time-dependent erodibility
To visualise the sensitivity of the modelled turbidity to the time-dependent erodibility function E β , we focussed on the main sampling site and computed * , the modelled turbidity time series with the value of E constrained to a constant value of * =1. The relative effect of erodibility on turbidity was then given by ( * ) . The results ( Figure 13) showed that the processes which we parameterised in the model as a seasonal variation in erodibility caused 50% attenuation of turbidity in mid-summer (greater attenuation for near-bed turbidity, less for near-surface). Removing an arbitrary value of half of the total available wave power at this site (averaged over the three years = 3.685 kW.m -1 ) would be equivalent to reducing the significant wave height to √0.5 = 0.71 of the natural state. Removing the same quantity of power by attenuating the tidal flow would represent only an 18% draw-down of the long term average current power, or a diminishing of the tidal speed to √0.82 3 = 0.936 of the natural state.

Sensitivity to changes in bed shear stress due to tide and wave energy extraction
Independently attenuating the significant wave height and the depth mean tidal current speed by these amounts and recomputed the bed shear stress, resulted in only small changes in predicted turbidity, and these were well within both the 95% prediction interval of the observed data in the natural system, and the 95% confidence intervals around the observed mean turbidity (Figure 14). Hence, we conclude that the impact of implemented decreases in either wave or current power on turbidity were unlikely to be detectable in the field with the available measurement capability.
The wave power resource at the main site was relatively small, so we also assessed the impact of removing a larger quantity of power purely by attenuating the tidal current speed (50% reduction in depth averaged current speed, corresponding to removal of 87.5% of the power). This resulted in modelled turbidity which was partially outwith the prediction interval of observed turbidity in the natural system, and completely outside the 95% confidence intervals around the observed mean ( Figure 14). Hence, we conclude that reductions in current speed of this magnitude would be detectable in the field.

Effects of tide and wave energy extraction on the 1% surface irradiance depth
Independently extracting 50% of the long-term average wave power or 18% of the average tidal current power, as outlined above, had an imperceptible effect on the underwater light environment at the main sampling site (long-term mean and s.d. of 1% irradiance depths:

Discussion
Simulation of suspended sediment concentrations is a notorious problem in coastal marine modelling, but is fundamental to many engineering applications and to improving understanding of marine ecology. The most difficult issues relate to a) the dynamics of vertical distributions of suspended particles in the water column under fluctuating flows, and b) the variability of seabed sediments with respect to their susceptibility to erosion. The former is exceptionally complicated given natural size distributions and properties of suspended particles, and flows arising from combinations of surface waves, tides and residual currents, but is at least in principle soluble from physical principles. However, there are no analytical principles to comprehensively predict the mobilisation of particulate material on the seabed. Commonly used formulations based on grain size and bed shear stress under constant flows appear to provide an explanation, but struggle to accommodate mixtures of grain sizes (El Ganaoui et al., 2004;Bartzke et al., 2013), fluctuating flows (Yu et al., 2011), bedform architecture (Soulsby and Whitehouse, 2005), sediment consolidation (McCave, 1984), and the variability induced by burrowing, sifting and habitat-modifying organisms (Amos et al., 1992). Hence, almost every model involves some massive assumptions and empirical parameterisations in order to render it applicable to a real-world situation. The virtue of involving detailed algorithms for some of some better understood aspects of the problem, as opposed to making simplifying assumptions when other aspects of the overall problem are only crudely implemented, is a difficult judgement. For example, the benefits of using turbulence closure schemes to model the vertical distribution of diffusivity as opposed to assuming generalised profiles may be minimal, when set against our approximate understanding of sediment erodibility.
Clearly, the transport of SPM is a three-dimension spatial problem, and if transport issues are the focus of attention then we have no alternative but to employ fully spatially resolved models, e.g. Gerritsen et al., 2000, review by Jones, 2002Mercier andDelhez, 2007, Danish Hydraulics Institute, 2013). However, all such models suffer from the fact that the spatial resolution of hydrodynamics is rarely if ever supported by corresponding spatial resolution of the parameters needed to characterise the seabed sediments. On the other hand, if, as in this study, the focus is on vertical concentrations of SPM then one-dimensional (vertical) models may be more appropriate, allowing detailed attention to hydrodynamics, the role of seabed sediment properties (Clark and Elliott, 1998;Dobrynin, 2009;Ramakrishnan et al. 2013). We used this approach here to develop a highly simplified parameter-sparse model capable of numerical optimisation, so as to capitalise on the substantial data resource available from the Stonehaven sites.
Our model incorporated a number of gross simplifications in pursuit of a parameter-sparse scheme. In particular, we characterised the seabed sediment only in terms of mud content (grains < 0.063 mm diameter) since we were primarily interested in SPM in the illuminated upper layers of the water column which would be expected to comprise only fine grained material. We also caricatured the well-known phase-lag between fluctuations in current velocities and SPM concentrations in the water column (Yu et al., 2011) by means of a timelagged average of the shear stress, and assumed a linear vertical profile of diffusivity. Nevertheless, our results show that the method was extremely successful in representing both the vertical distributions and dynamics of SPM concentrations at the time scales of the spring-neap tidal cycle and storm event. Fluctuations in SPM at sub-tidal (flood-ebb) time scales were also reproduced by the model, though probably with reduced amplitude compared to the observed data. As an independent credibility check, the fitted sinking rate of particles was consistent with the measured turbidity being due to grains in the 0.02 -0.03 mm size range which seems entirely reasonable. The model did not capture the extremes of high turbidity at the sampling sites, but this is to be expected since the current data used to derive the bed shear stress included only the tidal constituents, not any wind driven or surgedriven flows which would be expected to enhance the shear especially during storm events leading to higher extremes of turbidity. The most surprising aspect of the model parameter fitting was the implied magnitude of the seasonal effect which we modelled as a sediment consolidation process. The seasonal cycle of turbidity was very obvious in the observed data, and circumstantially one might assume that this could be explained by seasonality in sediment suspension due to wave action. However, the model shows that seasonal variation in bed shear stress is in no way sufficient to account for the seasonality in turbidity. Possible alternative explanations are seasonal stratification of the water column and modification of the vertical diffusivity profile, lateral transport related to seasonally varying inputs from river discharge, or seasonally varying seabed sediment consolidation. Summer thermal stratification of the water column in the study area was slight and highly transient, so this seems an implausible explanation for the seasonal cycle of turbidity. If we take salinity as an index of lateral transport from the nearshore environment to offshore, we can see that salinity in the study area also followed a distinct seasonal cycle but around 4 months out of phase with the seasonality in turbidity. Minimum salinity off the east coast of Scotland occurs in April (Bresnan et al., 2008), coinciding with the peak in discharge from the main rivers (Dee and Don) whose catchments include the mountainous regions of the Scottish Highlands. In contrast, the fitted seasonal cycle of consolidation was more or less symmetrical around mid-summer/mid-winter, consistent with observed data on the algal pigment content and microbial respiration rates of sediment cores (Serpetti, 2012). Hence, we conclude that our inclusion in the model of seasonal consolidation was a credible explanation for the seasonality of turbidity, rather than any advective effect.
Biological consolidation of sediments may arise through a variety of processes. Secretion of sticky organic molecules by microbes Gust, 1987, Lubarsky et al., 2010), benthic algae and microbes clogging the pore spaces and binding grains together (Nowell et al., 1981;Sutherland et al., 1998;Austen et al., 1999;Paterson and Black, 1999), and forming mats on the sediment surface all lead to inhibition of sediment erosion (Fonseca, 1989;Paterson, 1989;Oppenheim and Paterson, 1990). Living algal mats are most prevalent in shallow waters since the micro-organisms concerned require light to photosynthesise. However, other biological processes may have the opposite effect on sediment erodibility due to destabilisation of the sediment structure. These include bioturbation by burrowing and sediment ingesting macrofauna and meiofauna which reprocesses sediment into faecal granules (Montague, 1986;Amos et al., 1992;Rowden et al., 1998;Lumborg et al., 2006).
A key question concerns the transferability of parameters for our model to other locations where the depth, tide, wave and sediment characteristics may be very different from those off Stonehaven. The assumption of a linear profile of vertical diffusivity and the Rouse formulation of declining turbidity with altitude in our model may be inappropriate in extreme high flow situations, but in these conditions the seabed sediment mud content should be negligible so the mass of fine-grained material in suspension will be small. Any suspended material is likely to be coarse grained with a fast sinking rate, and contribute relatively little to turbidity at high altitudes above the seabed. However, we should certainly expect some regional specificity of some of the seabed sediment parameters, in particular the sensitivity parameters for mud content (ε) and erodibility (β). Since we do not model horizontal advection of SPM we cannot explicitly take account of the exchange of SPM between adjacent sites of different seabed mud content, which must occur in reality. However, these effects are implicitly included in the term S ε . We might expect ε to approach 1 for sites within a large homogeneous area of seabed. In contrast, the Stonehaven study area is a complex network of sediment patches of different grain size composition at length scales equivalent to the tidal excursion (Serpetti et al., 2011), so we would expect turbidity over patches of coarse sediment to be higher than in a homogeneous system, with values of ε <<1.
Similarly, the extent of biologically induced consolidation and erodibility is highly likely to be region-specific. The phenomenon is well known and extensively studied in tidal mud-flats and shallow estuaries where the sediments are predominantly fine cohesive muds and the effects of biological activity are very obvious (Le Hir and Karlinkow, 1992;Austen et al., 1999;Paterson et al., 2000;Widdows et al., 2000;Andersen, 2001). In fact, seasonal variation in erodibility mediated by biological activity may be the dominant factor controlling water turbidity in shallow tidal regions such as the Wadden Sea (Lumborg et al., 2006;Borsje et al,. 2008;De Vires and Borsje, 2008,). In contrast, models of sediment suspension and transport in deeper open shelf systems generally assume that spatial, and especially temporal, variability in erodibility due to biological consolidation can be disregarded (e.g. Pohlmann and Puls, 1994;Ribbe and Holloway, 2001;Kuhrts et al., 2004;Pandoe and Edge, 2004;van den Eynde, 2004). However, recent research shows that this may not be the case (Stevens et al., 2007;Briggs et al., 2015). For example, Dobrynin (2009) found that a model of suspended sediment concentrations in the southern North Sea was unable to explain the distribution of surface concentrations derived from satellite remote sensing without resorting to alternative summer and winter parameterisations of erodibility.
Taking the issues of sediment heterogeneity and erodibility together, it is likely that our model is still not general enough to be parameterised in one region and directly applied in another. The limitation on re-calibration for a new regional application is likely to be the availability of data on vertical profiles of turbidity which are rarely measured as part of routine oceanographic surveys. There are several ways in which we could address this. First would be to make use of satellite remote sensing data on SMP concentrations (e.g. Dobrynin, 2009;Rivier et al., 2012;Sabatino et al., 2015). These do not immediately provide vertical profiles of SMP concentrations, but they do offer near-surface horizontal distributions across contrasting seabed environments and water depths. Secondly, we could explore measurements which might be correlates of erodibility and spatial heterogeneity of sediments, and incorporate these into the model. For example, the algal pigment content of sediments has been investigated as potential indicators of biologically-mediated consolidation (Riethmuller et al., 2000), but so far has not shown general applicability.
The primary motivation for our study was to scope the regional impact of tidal and wave energy extraction on water column turbidity. For the main sampling site in our study area, we show that removing an annual average of around 4 kW.m -1 by either attenuating the significant wave height to 70% of the natural state, or the tidal current speed to approximately 90% of the natural state, produces reductions in turbidity and light penetration into the sea that would likely not be detectable given the variability in measurement. These proportional attenuations of wave height and tidal current speed could be regarded as realistic expectations for the 1-20 km scale footprint of large scale energy extraction arrays (e.g. Wu et al., 2015).
Our results show that much larger attenuations of current speed, in the order of 50%, are required to have effects on turbidity and underwater light climate which might be detectable at least in our sampling region. Flow changes of these magnitudes may be expected in the immediate vicinity of extraction devices such as turbines, but then other wake and small scale turbulence effects not represented in our model will probably dominate. However, the problem with these scoping estimates is that they assume the seabed sediment characteristics are fixed and independent of changes in the flow regime. Closing this feedback connection between alteration of the seabed sediment landscape due to changes in hydrodynamics, and the supply of material for suspension into the water column, is a significant modelling challenge. For example, simulations of the hydrodynamic implications of introducing turbine arrays in the Pentland Firth (UK) and the Bay of Fundy (USA) indicate significant accelerations and decelerations of the flow regime over a large area, leading to changes in bed-load transport and alteration of the sediment distribution (Fairley et al., 2015;Martin-Short et al., 2015;Wu et al., 2015). Our future aim is to use such predictions of the altered sediment landscape, together with the altered hydrodynamics, to predict the 1-dimensional (vertical) distributions of turbidity and light environment at a network of spatial locations.
available at DOI:10.15129/6f7f7a5f-193f-475c-9568-1c0c18b015a5, and turbidity data at DOI: 10.15129/ea4d4dbf-9807-4063-a4e2-716a1a6b7f7b.         Wu and Wang, 2006), and the sinking rate estimated by the model fitting process (horizontal blue lines and grey shading). The grey shading indicates ±1 s.e. around the fitted sinking rate value. The vertical dashed line to the right indicates the upper grain size threshold for mud. The fitted range of sinking rates implies a suspended grain size of 0.02 -0.03 mm.

Fig. 7.
Average vertical profiles of modelled and measured turbidity at each of the sampling locations. Heavy red lines represents the median of the observations at each 0.5m interval of altitude above the seabed at each site, whilst the red dashed lines span the 5 th -95 th centiles. Heavy black line at each site represents the median of the modelled turbidity for the subset of model output times corresponding to the times of the observations, while the grey shading spans the 5 th -95 th centiles.     Model sensitivity to seabed mud content. The plotted line represents turbidity relative to that at the main sampling site as the mud content of seabed sediments is varied between a small value >0, and 1. The symbols show where the other sampling sites lie on this relationship. Fig. 13. Model sensitivity to seasonal erodibility. The plotted lines shows the turbidity at two depth horizons at the main sampling site during 2009, relative to a scenario in which the erodibility term in the model was held constant and equal to 1.

Fig. 14.
Model sensitivity to extraction of wave and tidal energy. X-axis refers to the modelled turbidity in the natural system at the main site. Y-axis refers to either the observations of turbidity in the natural system, or modelled scenarios of energy extraction. The diagonal heavy dashed line corresponds to a 1:1 relationship which would represent a perfectly fitting model of the natural system, or extraction scenarios causing no change in turbidity. Large symbols show the observed data from the natural system (open symbols 5 m depth, filled symbols 5 m altitude). The thin dashed lines correspond to the 95% confidence interval around a regression of modelled vs observed natural state data (inner pair of lines), and the 95% prediction interval of observed turbidity given modelled values (outer pair of lines). Very small symbols correspond to modelled turbidity at depth horizons given energy extraction scenarios: black: 5 m altitude, 3.7 kW.m -1 wave extraction; blue: 5 m depth, 3.7 kW.m -1 wave extraction; red: 5 m altitude, 3.7 kW.m -1 tidal extraction; green: 5 m depth, 3.7 kW.m -1 tidal extraction; grey: 5 m depth and 5m altitude, 50% attenuation of tidal current speed.