MERAMOD : predicting the deposition and benthic impact of aquaculture in the eastern Mediterranean Sea

A model composed of coupled particle tracking and benthic response modules was developed for predicting waste solids flux and benthic impacts of gilthead sea bream Sparus aurata L. and sea bass Dicentrarchus labrax L. aquaculture. The model was tested at 6 sites with different hydrodynamics, bathymetries and biomasses in the Aegean and Ionian seas, eastern Mediterranean Sea, and sediment trap flux and benthic impact indicators were observed. Seven sediment trap validation studies were conducted that varied in design with traps deployed either on the seabed, attached to nets or suspended in the water column. Model predictions of flux to traps spaced 5 m apart up to 50 m from the cages over a 13 d period were significant (R2 = 0.61, n = 57, p ≤ 0.05). However, the model could not predict adequately the flux to traps spaced 2 m apart in the high-flux zone underneath cages where variability between trap observations was high. In this high-flux zone underneath cages, the averaged model flux predictions resulted in a performance of ±49%. Statistically significant relationships were established at 4 sites between modelled flux and either benthic fauna impact indicator species (S), abundance (A), A/S ratio, Shannon-Wiener index or biomass fractionation index (BFI), (R2 = 0.82, 0.60, 0.57, 0.67 and 0.48, respectively; n = 24, p ≤ 0.05). Two other sites, which did not exhibit an abundance peak in enriched zones, did not fit these relationships. Using relative abundance of taxonomic groups, a modelled flux of 4.1 g m−2 d−1 was determined to be a useful boundary; on either side of this boundary, clear trends occurred in pollutant tolerant and intolerant species.


INTRODUCTION
In the 20 yr from 1989 to 2009, Mediterranean aquaculture of marine finfishes rose from 3384 to 222 149 t; 95 677 t were produced by Greece and77 196 t by Turkey in 2009 (FAO 2011).The benthic impacts of the Mediterranean fish farming industry are still less well described and understood than those at higher latitudes, where salmon is the pre-dominant culture species group, but considerable progress has been made in recent years in understanding and quantifying the relationships between particulate waste inputs and benthic response in this environment (e.g.Holmer et al. 2008, Apostolaki et al. 2009, Papageorgiou et al. 2009, 2010, Sanz-Lázaro et al. 2011).Regulation and management of seabed impacts is a key component of regulating finfish impacts in all salmon farming countries (Wilson et al. 2009) but is inconsistently applied in Mediterranean aquaculture (Telfer et al. 2009).
Waste deposition models that predict flux of waste material at the seabed are relatively common (Hevia et al. 1996, Silvert & Sowles 1996, Silvert & Cromey 2001, Cromey & Black 2005, Jusup et al. 2009, Symonds 2011), but very few models relate flux to a benthic effect.Aquaculture impact models that link the physics (flux) to a benthic effect (biological) and/or include a biochemical component are useful for planning and regulation of aquaculture sites, as measurement of benthic community structure is a well-established method of assessing the effect of mariculture discharges (Brown et al. 1987, Findlay & Watling 1997, Pearson & Black 2001, Brooks & Mahnken 2003, Pereira et al. 2004, Borja et al. 2009).In the present paper, we used a large body of benthic impact and sediment trap data collected from several fish farms in the eastern Mediterranean Sea to develop a computer model that will allow prediction of benthic impacts and thereby aid good environmental management and regulation.
DEPOMOD, a model previously developed to predict benthic impact from salmon farms (Cromey et al. 2002), was generally found to underpredict the deposition and benthic impact of aquaculture when tested at several fish farms in the Mediterranean (authors' unpubl.data).Wild fish aggregate in large numbers around Mediterranean fish farms (Dempster et al. 2002, 2004, 2005, Dempster 2005) to scavenge for waste particulate material.However, this process was not incorporated in DEPOMOD applied at salmon farms even though wild fish are now known to aggregate around these northern temperate fish farms (Dempster et al. 2009).Magill et al. (2006) determined the difference between the settling rates of faecal particles from gilthead sea bream Sparus aurata L. and sea bass Dicentrarchus labrax L. and demonstrated the importance of using these rates in modelling studies; fish farms in the eastern Mediterranean Sea typically have both sea bass and gilthead sea bream from different cohorts on the same farm, where feed input may differ by an order of magnitude in adjacent cages.
The objectives of the present study were to develop and validate a model capable of predicting the benthic impact of deposition from gilthead sea bream and sea bass aquaculture in the eastern Mediterranean Sea.It was necessary to develop the capability to specify faecal settling rates on a cage by cage basis in the new model (henceforth MERAMOD), as well as to specify a distribution of faecal settling rates representative of the coarse, fast-settling faecal par-ticles from sea bass.Thus, (1) new processes were developed in the model specific to the cultured species, and (2) species-specific input data for sea bass and gilthead sea bream were determined.To calibrate and validate MERAMOD an extensive programme of sediment trap studies and spatially referenced benthic sampling investigations was carried out in the Mediterranean Sea.The aim of the benthic response model was to establish relationships between modelled particle flux to the seabed (total dry solids [TDS] as g solids m −2 yr −1 ) and benthic community descriptors so that MERAMOD can be used for planning and monitoring Mediterranean aquaculture.

Time and study sites
Site notation is given as aMDb, where a is the sampling time (1 = July 2001, 2 = March 2002, 3 = October 2002), b is the site reference number (1 to 6, 8), and 2MD1A and B represent experiments conducted on the same site and month on different days.This anonymous site reference system is used as the data were collected under a confidentiality agreement.
To compare observed and predicted particle flux, 2 types of sediment trap studies were undertaken: short period deployments (24 h) in high-flux zones (i.e.directly underneath cage groups) and one deployment over a 13 d period incorporating low-flux zones sufficiently far from the cages to receive little or no farm input (3MD8 study site).Traps were deployed on the seabed, attached to the nets and suspended in the water column (6 experiments with 24 h deployments and 1 experiment lasting 13 d; Table 1, sites prefixed 2 and 3 only).The benthic response model was validated using benthic faunal data collected along transects from 6 sites (Table 1, sites prefixed 1).

General model set-up
Grid cell resolutions of 10 m and 1 m were used for benthic response model validation and sediment trap validation, respectively, which reflected the spacing of the stations along the transects (13 d study at 3MD8) or replicate traps (24 h studies at MD1, MD3 and MD5).Particle steps were evaluated every 6 s to give a vertical step length of 6 cm for a particle settling at 1 cm s −1 .Flux predictions at specific points in the model grid (e.g.trap positions) were determined by using nearest neighbour interpolation.In sediment trap studies, particle starting positions in the cage volume were dependent on the degree of fouling and the pattern of flux determined from traps attached directly to the nets.For the longer term simulations of 1 mo, random starting positions of particles were assigned inside the cage volume.In all simulations, total particle numbers were optimised (> 5 × 10 5 ).
Time steps of 5 min (equivalent to the frequency of measured current data) were used for trap validation studies, and hourly averaged data was used for the longer benthic response model runs.A sensitivity analysis demonstrated the need to use all water current data without averaging for short-term model runs, whereas longer runs of a month or more were represented satisfactorily with hourly averaged data.In all runs, faecal settling data for sea bass of 0.4 (6%), 1.4 (9%), 2.5 (20%), 3.6 (38%) and 4.6 cm s −1 (27%) and for gilthead sea bream of 0.4 (24%), 1.5 (45%), 2.5 (18%), 3 cm s −1 (13%) were used (Magill et al. 2006).Where species in individual cages were unknown, as was the case when only monthly summaries were supplied by the farmer, a combined distribution of settling rates for both species was used: 0.4 (8%), 1.4 (14%), 2.5 (20%), 3.5 (35%) and 4.6 cm s −1 (23%).
A general figure of feed wastage (F wasted ) of 5% was used for all sites, and feed digestibility (F dig ) and water content (F w ) were set at 85 and 9%, respectively.Site-specific values for these variables were not determined, and default values were applied consistently across sites.Accurate determination of long-term feed wastage from marine cages is technically very difficult and is likely to vary from site to site and over time.We chose 5% for F wasted as this value has given satisfactory model predictions in other contexts (Cromey et al. 2009).The feed digestibility of 85% is taken from feed manufacturers' technical data but likely varies between brands and over time, and the 9% water content was estimated by drying and weighing feed samples.
Feed input was obtained from the farmer at the time of each survey.For 24 h trap deployments, detailed real-time feed input data were obtained for the study cage and surrounding cages.For the benthic response model, average feed input data were used for the 3 mo before benthic sampling began.These monthly averaged data were typically for the whole farm, but occasionally cage group data were obtained.Where whole farm summaries were used, detailed husbandry data obtained for sediment trap studies (March 2002, October 2002) were used to scale farm summaries so that the feed input in cages above benthic sampling transects were represented more accurately.Empty cages were accounted for where these were re corded.
S B a Combined distribution for 6.5 and 5.0 mm pellets; b combined distribution for 2.5 and 5 mm pellets Table 1.Input data for benthic response model (sites with prefix 1 for July 2001) and sediment trap studies (with prefixes 2 and 3 for March 2002 andOctober 2002, respectively).Species are mixed, sea bass (Ba) or gilthead sea bream (Br).Fish weight given if known.Faecal particle release times are I: continuous; II: 09:00-10:00 h (40%), 60% remainder; III: 10:00-11:00 h (35%), 15:00-16:00 h (35%), 30% remainder; IV: 08:00-16:00 h (70%), 30% remainder; V: 09:00-18:00 h (70%), 30% remainder.Deployment methods: benthic response study -seabed deployments of 0, 5, 10, 25 and 50 m from the cage edge in addition to 2 distant reference stations (Bn); sediment trap studies -traps attached to underside of net (NET), on the seabed (SB) or suspended in the water column (WC).For deployment times, see Table 2.For trap studies, information is only shown for the study cage (see 'Materials and methods' for explanation of trap deployments).nd: no data available Feed settling data (described in 'Feed settling velocity' below) were specified depending on the pellet diameter and brand, and when this information was not available, a generalised relationship for feed settling velocity was used.In sediment trap studies, diurnal patterns of release of feed and faecal particles were matched with husbandry activities.This varied from 1 or 2 feeding events per day to continuous feeding in daylight hours.However, for the longer benthic response model predictions, waste faecal and feed material was released continuously as the specific timing of feeding and defecation events was not sensitive at this time scale.

Bathymetry and hydrography
Bathymetry and cage positions were measured at each site with the RoxAnn™ system interfaced with differential GPS (dGPS) (MacDougall & Black 1999).Positional data were corrected to Universal Transverse Mercator (UTM) coordinates (zones 26S, 27S), but depth corrections were not necessary as all sites were micro-tidal (tidal range, < 30 cm).Data were interpolated by using the kriging algorithm (slope = 1, anisotropy ratio = 1 and angle = 0, output grid spacing = 10 m) and incorporated into MERAMOD, where the modelled grid sizes were sufficient to include the whole initial deposition footprint (grid areas range, 0.17 to 0.26 km 2 ).
To measure water current velocities, electromagnetic current meters (S4, InterOcean Systems) or an Acoustic Doppler Current Profiler (ADCP) (RDI 600 kHz Workhorse, Teledyne RD Instruments) were deployed during trap deployments (Table 2).Three S4 instruments were deployed at the surface, net bottom and nearbed hanging from the cages and had a sampling interval of 5 min (2MD1A, 2MD1B, 2MD5 and 3MD8), or in the case of the ADCP, deployed near the surface looking down (3MD1, 3MD5 and 3MD3; bin size = 1 m).
To obtain longer time series for validation of the benthic response model, hydrographic data were measured at 4 sites (MD1, MD3, MD5 and MD6) around the time of benthic sampling (July 2001), and 2 sites had data available for other periods (MD2, December 1998;MD4, July 1996).MD2 was the only site where data were not available for the summer season (the season of benthic sampling).S4 or rotary current meters (RCM7, Aanderraa Instruments) were deployed on a U-shaped mooring within 100 m from the cages at the surface and near-bed (except at the shallow site 1MD4, where the deployment was at mid-depth).
A meteorological station (Weather Wizard III, Davis Instruments) was deployed concurrently at nearly all sites either on cages or on land, depending on security, with the anemometer cup situated 3 m above the land or sea surface.A sampling interval of 10 min was used for 24 h sediment trap studies and 20 min for longer deployments.These data were used to analyse the influence of winds on surface currents and were not used directly in MERAMOD.
Drifting buoys (drifters) deployed at the 6 sites used for validation of the benthic response model (Sites MD1−MD6) were used to measure horizontal dispersion coefficients for use in modelling (Yanagi et al. 1982).At each site, a base station was located so that a line of sight could be achieved for all possible drifter positions.The drifters consisted of a GPS aerial and receiver and an ultra high frequency (UHF) transceiver, powered by a 12 V lead acid battery.Dimensions of each buoy were 17 cm in diameter by 29 cm with a mast height of 37 cm, and each buoy had a foam collar to provide buoyancy.The exact position of the base station was determined by averaging GPS measurements for 30 min so that subsequent GPS readings received by the base station could be corrected for inherent errors.Corrections were then applied to determine the dGPS position of each buoy (accuracy: 57% data ± 1 m, 99% data ± 4 m).Six buoys were deployed together, each transmitting its position every 30 s.The drogue sock length was 2 m, and the sock middle was situated 6 m beneath the surface.Deployment length varied from 1 to 7 h.Dispersion coefficients were calculated according to Yanagi et al. (1982, see Appen dix 1).Dispersion varies continuously over time and space, and drifters were retrieved when they were advected more than 1 km from the site or snagged on mooring equipment.Dispersion coefficients, k x and k y , calculated from the longest deployment at each site, were used in the modelling.Vertical dispersion coefficients were not measured and assumed to be 0.001 m 2 s −1 (Gillibrand & Turrell 1997).

Feed settling velocity
Settling velocities of commercially available sea bass and sea bream extruded feeds were determined by timing the distance to settle over a known distance in a settling column.Settling velocities were measured over a distance of (1) 42.5 cm in a glass column (internal diameter [Φ] = 8.5 cm, length = 70 cm, temperature = 26.5°C,salinity = 37.6) and (2) 1.00 m in a 4 mm thick Perspex column (internal Φ = 15 cm, length = 2 m, temperature = 26.2°C,salinity = 37.1).
In these experiments 11 pellet types were used (Φ = 0.9 to 6 mm).The 2 different length columns were necessary as experiments were undertaken in different geographical locations, but this did not affect the results as theory suggests that 6 mm pellets settling in seawater reach terminal velocity in less than 10 cm (Cromey et al. 2009).

Model validation: particle tracking
Modelling was undertaken with detailed information for each cage including species, feed inputs for the time of survey and faecal settling velocities.Sensitivity analysis showed that sedimentation under-neath a cage is influenced by surrounding cages as well as the cage above, so detailed information on the surrounding cages was essential.The wild fish module settings were varied according to observations at each site and, in the case of 3MD8, it was used to test the effect on model performance by examining regression line parameters (optimise for gradient closest to 1 and intercept to 0, described in 'Mass balance checks and statistical testing of model performance' below).
Detailed information on the degree of cage fouling and feeding technique was also obtained, and, at some sites, a description of where food pellets entered the cage at the surface.In some of the runs, this level of detail was included in the model.Predic- tions of deposition (g m −2 ) were obtained over the modelled period and then standardised by scaling up to obtain flux (g m −2 yr −1 ).PVC sediment traps met the well-accepted criterion of length:diameter ratio of at least 5:1 (Wassmann et al. 1991) and, in most cases, traps were deployed by a diver.At site 3MD8, divers deployed sediment traps (internal Φ = 9.6 cm) for a period of 13 d directly underneath the study cage centre and then on 8 radial transects at distances of 5,10,15,20,25,35 and 50 m from the centre (a total of 57 traps).The contents were analysed for ash-free dry weight (AFDW).For all other trap studies, traps (internal Φ = 15.4 cm) were deployed for 24 h in 2 m intervals either directly attached to the net bottom (denoted NET), suspended in the water column by rope attached to the cages (denoted WC), on the seabed (denoted SB) or a combination of these (Table 1).WC trap deployments consisted of 12 traps on a line deployed across the width of the cage from above approximately 2 to 5 m below the base of the net, whereas SB and NET trap deployments consisted of 6 traps each.In these short-term experiments, a study cage was selected that contained large fish to maximise the signal from the farm.By using different deployment depths, the robustness of flux predictions at different water column depths was thereby tested.
At 3MD8, background sedimentation increased in importance compared with the farm flux for more distant stations and was determined from 3 sediment traps deployed (>100 m) outside the estimated footprint zone of deposition (predicted by pre-survey modelling).During 24 h experiments at the other sites in the high-flux zone (>10 000 g m −2 yr −1 ), background sedimentation was assumed to be negligible compared with the very high-flux zones beneath the cages, so reference traps were not necessary.
Traps were retrieved, contents allowed to settle until clear (~15 min) and excess water siphoned off while care was taken to not disturb the contents.To obtain TDS, contents were then filtered through preweighed coarse glass fibre filters with a vacuum apparatus (Whatman GF/A, pore size 1.6 μm, Φ = 15 cm) and then dried at 60°C overnight.For the longer deployment at 3MD8 only, we ob tained AFDW (AFDW = organic matter [OM]) by combusting the sample in a muffle furnace at 520°C overnight and subtracting the ash weight from the TDS weight.In both determinations, samples were allowed to cool in a desiccator before they were weighed to an accuracy of 0.001 g.By using the trap diameter, observations were scaled to obtain standard units of g m −2 yr −1 to compare with model predictions.

TDS flux and benthic impact indicators
Flux predictions were obtained at each sampling station for which observed benthic community descriptors were available.A logarithmic scale (modelled flux + 1) was used on the x-axis so that undercage (TDS flux, >10 000 g m −2 yr −1 ) and reference stations or outlying stations could be included on the same scale.A logarithmic scale was used for abundance and biomass descriptors on the y-axis by convention.From diver observations, wild fish feeding on waste pellets was an important process at all sites.
Benthic fauna sampling By using current meter data, a transect was established at Sites MD1 to MD6 along the prevailing direction of the water currents and close to high biomass cages.Samples were taken at a distance of 0 m (under the cages) and at 5, 10, 25 and 50 m from the edge of the cages downstream.Reference stations were chosen approximately 1 km upstream from the cages at similar depth and sediment type.Samples at 0, 5 and 10 m from the cages were taken by SCUBA divers with Plexiglas cores (internal Φ = 9.5 cm, penetration depth = 15 cm) while at 25 and 50 m and the reference station, samples were taken by means of a Smith-McIntyre grab (sampling area = 0.1 m 2 ).
At each station, 5 replicate samples for macrofaunal analysis were taken.All macrofaunal samples were sieved in situ through a 500 μm mesh, and the retained sediment containing macrofaunal organisms was preserved in 10% buffered formalin.In the laboratory, macrofauna were sorted, enumerated and identified to major taxa, i.e.Polychaeta, Crustacea, Echinodermata, Mollusca, Sipuncula and miscellaneous (including Anthozoa and Nemertina) under a stereomicroscope.The wet mass from each taxon was measured to the nearest 0.1 mg after blotting the sample dry on filter paper, and, after biomass measurements, all taxa were further identified to species.Indicators extracted from the benthic data set included macrofaunal abundance, biomass, Shannon-Wiener diversity index (H ') and the biomass fractionation index (BFI).The BFI is defined by Lampadariou et al. (2008) as the macrofaunal biomass retained on a 0.5 mm sieve only, divided by the macrofaunal biomass retained on 1.0 mm plus 0.5 mm sieves, when macrofauna are sieved sequen-tially over 1 mm and 0.5 mm sieves.Low BFI values indicate a low proportion of opportunistic (typically small) species.

Indicator species and families
The indicator families and species were derived by applying a similarity percentages (SIMPER) analysis routine (Clarke 1993) of the PRIMER software package on the macrofaunal data set of all the sites and stations.To ensure that the indicators chosen were representative of the eastern Mediterranean Sea in general, all taxa that did not occur at all the sites sampled were excluded.To select indicators typical of disturbed conditions near the cages and of relatively undisturbed conditions distant from the cages, all taxa that occurred at all stations at each site were excluded as well.
For each station, the abundance, relative abundance and presence or absence of indicator species and families were calculated.Indicator species of polluted conditions used were Apseudes robustus, Capitella capitata, Caulleriella oculata and Maldane sarsi.Indicator species of unpolluted conditions used were Cirrophorus branchiatus, Cossura coasta, Levinsenia gracilis, Magelona alleni and Monticellina heterochaeta.Indicator families of polluted conditions were Ampelscidae, Lucinidae, Maldanidae, Nemer tina and Oedicerotidae.Indicator families of unpolluted conditions were Ampharetidae, Onuphidae, Opheliidae and Semelidae.Relationships between predicted flux and the relative abundance of indicator species and families were then established.

Sediment chemistry
Measurements of redox potential, sediment particle size, total organic carbon (TOC) and nitrogen and plant pigments were performed by standard methods as previously described in Lampadariou et al. (2008).As in that work, we report redox at 4 cm depth after Pearson & Stanley (1979).

Mass balance checks and statistical testing of model performance
A spreadsheet calculation was used to check that the total mass accumulated in the grid equalled the solids discharged from the cages (see Appendix 1).
For the sediment trap studies, calculations in a spreadsheet allowed some assessment of the quality of the observed data, particularly for the NET traps, which may be subject to receiving biodeposits from the farm such as encrusting biota, which are not part of the waste feed and faecal components.On one occasion non-fish farm material was observed at the time of filtering and these data were rejected (see Appendix 1).As NET traps were often deployed across the width of the net bottom, these basic calculations also revealed whether the waste material was exiting the net evenly across the whole area of the net bottom or in some other manner.This important detail was incorporated into some modelling scenarios and was caused by uneven distribution of feed and faeces in the nets resulting in some traps re ceiving a greater flux than adjacent traps.
Models were assessed with the linear regression of Y i = β 0 + β 1 X i for observed (Y i ) and predicted (X i ) values, where β 0 and β 1 are the intercept and gradient, respectively.Plotting predicted values on the x-axis is widely accepted, as model predictions are more under the control of the modeller and may be deemed to contain less error than observations from an irregular sampling regime (Sokal & Rohlf 1995).A Student's t-test was undertaken for R 2 , β 0 , β 1 (α = 0.05, df = n − 2) to test if R 2 , β 0 and β 1 were significantly different from 1, 0 and 1, respectively (Mesple et al. 1996, Jusup et al. 2009).Following the termino logy given in Jusup et al. (2009) (after Tett 2007), models were classified as follows: Excellent, the slope and intercept are not significantly different from 1 and 0, respectively, and R 2 is sig nificantly greater than 0; Good, either the slope is significantly different from 1 or the intercept is sig nificantly different from 0, and R 2 is significantly greater than 0; Fair, the slope is significantly different from 1 and the intercept is significantly different from 0, and R 2 is significantly greater than 0; Poor, the variance ex plained (R 2 ) is not significantly different from 0.
In the benthic response model, curves were fitted on the basis of minimising for each comparison a measure of the goodness of fit: Σ(O − E) 2 /E where O is the observed value of indicator and E is the modelled value of indicator.In addition, R 2 was also maximised.
Model performance was also expressed by averaging the percentage difference between predicted and observed flux.For example, for a predicted flux of 1000 g m −2 yr −1 with a specified performance of ± 20%, observations would be expected to fall in the range of 800 to 1200 g m −2 yr −1 .

Bathymetry and hydrography
The study site transects ranged between 11 and 31 m depths, and the mean surface current speeds varied between 1.8 and 4.1 cm s −1 as determined from the longer current records.Dispersion coefficients varied between 0.003 and 0.42 m 2 s −1 and did not correlate with surface current speed.Model predictions are somewhat sensitive to horizontal dispersion coefficients: using Eq.6 (Appendix 1) and a 60 s interval (δt) a k x of 0.1 m 2 s −1 gives a random step of 3.5 m compared with 7.0 m when k x is 0.42 m 2 s −1 .Use of survey data for the longest runs was thus preferable to using the default value of 0.1 m 2 s −1 used in DEPOMOD.In addition, by taking measurements, we were able to specify k x and k y according to the site observations, rather than using an identical k x and k y as with a default setting.
By considering surface current for all measurements and dispersion coefficients, the following ranking of sites was determined (most dispersive first): MD8 > MD6 > MD4 > MD5 > MD2 > MD1 > MD3 (although none of the sites can be regarded as dispersive in absolute terms), and ranking was also determined in terms of depth (deepest first): MD5 > MD3 = MD6 > MD2 > MD1 = MD8 > MD4.The influence of wind on current speed and direction was unpredictable at many of the topographically constrained sites.For example, at 1MD5 surface currents were to the NE, but winds were from the NNW at first and then from the NNE (record length, 45 d; 10 min sampling interval), i.e. the surface current flowed against the prevailing wind direction.This indicates that representative current time series are essential for accurate modelling at these micro-tidal sites.

Feed settling experiments
Mean settling velocity (vs) increased with mean pellet diameter (Φ) up to 4 mm (Fig. 1), although vs did not increase further for 5 and 6 mm pellets (n = 455).Thus, as pellet density was not constant with diameter, use of these measured data in modelling is more accurate than applying Stokes' law to different particle diameters.The mean and SD of the whole data set were 10.9 and 3.8 cm s −1 , respectively, and these values were used in model scenarios where pellet diameter was unknown.

Model validation: sediment traps
Comparisons between observed and predicted AFDW at MD8 resulted in an 'Excellent' model fit (sensu Jusup et al. 2009, see 'Materials and methods: Mass balance checks and statistical testing of model performance' above; regression line: Y = 1.04X + 82 g m −2 yr −1 , R 2 = 0.61, p = 0.05, n = 57; Fig. 2d) after subtracting background sedimentation from observations for each sampling station and when 50% of waste feed pellets were removed by wild fish feeding.When considering data from all 57 seabed traps deployed at 3MD8, the mean difference between observed and predicted values was ± 44% of the predicted value (Table 3).However, the model performed better in the mid-and high-flux zones (± 29% and ± 35%, respectively) than in the low-flux zone (±111%, Table 3).The statistics of the regression lines shown in Fig. 2 are shown in Table 4.
In the 24 h experiments, observed and predicted fluxes were averaged for all traps in each under-cage deployment (Fig. 3), as the high variability observed between traps spaced 2 m apart was not reproduced by the model.Generally, both observed and predicted fluxes were higher in the NET traps compared with SB traps, with the exception of 2MD1A, where adjacent cages with higher feed input may have influenced flux to the SB traps.Model performance was higher for SB traps (classed as 'Good') compared with NET and WC traps (classed as 'Poor') (Tables 3  & 4), which was encouraging as model performance is important for predicting longer-term benthic impacts at the seabed.
In the experiments in the high-flux zone, observations and predictions exceeded 10 000 g m −2 yr −1 TDS at all stations, and maximum fluxes were greater Fig. 1.Mean (± SD) settling velocity of commercial extruded feed pellets used for sea bass and gilthead sea bream.Perla, Ichthys and Europa are brand names of fish feed used at the study sites than 100 000 g m −2 yr −1 (274 g m −2 d −1 ).Model performance was similar for the NET and WC traps at ± 39% (n = 18) and ± 36% (n = 36), respectively (models classed as 'Poor'), and better performance was evident at the SB traps at ± 66% (classed as 'Good') (Table 4).

Description of benthic impact
Alterations in the sediment chemistry occurred at sites in the proximity of the 6 fish farms investigated (MD1−MD6), although they varied in extent (Fig. 4) (see Lampadariou et al. 2008 for a preliminary assessment).Site 1MD4 was characterised by silty sediments with an average median particle diameter of 0.05 mm and a silt and clay content of > 60%, whereas Sites 1MD1, 1MD5 and 1MD6 were characterised by fine sand with an average median particle diameter of 0.30 mm and a silt and clay content of > 40% (Fig. 4a).Conversely, Sites 1MD2 and 1MD3 were characterised by coarse sandy sediments (median particle Φ = 0.60 mm, silt and clay < 7%).Redox values (Fig. 4b) at the silty seabed of Site 1MD4 were negative at all stations of the transect and increased above zero only at the reference station, whereas redox in the coarse sediments of Sites 1MD1 and 1MD2 remained positive even at 10 m distance from the cages.Sites 1MD3, 1MD5 and 1MD6, on the other hand, showed negative redox values up to 10 m distance but increased above zero at greater distances.Sites 1MD3, 1MD4 and 1MD5 showed increased values (up to ~5%) of sediment organic carbon content at the 0, 5 and 10 m stations indicating enriched conditions near the cages (Fig. 4c).In contrast, low carbon values were found at all stations at Sites 1MD1, 1MD2 and 1MD6.Macrofaunal abundance underneath the cages was relatively low and increased significantly at stations up to 10 m distance from the cages by a factor ranging from 5 to 25 for Sites 1MD1, 1MD2, 1MD3 and 1MD6 (Fig. 5a).However, at distances greater than 10 m for these 4 sites, macrofaunal abundance decreased rapidly and was comparable with values observed at the reference stations.On the contrary, abundances up to 50 m distance from the cages at the other sites (1MD4, 1MD5), were relatively low and similar to the reference stations.At Sites 1MD1 and 1MD5, biomass was impoverished up to 10 m distance from the cages, whereas it increased significantly and above the levels of the reference stations at the 25 and 50 m stations (Fig. 5b).For the remaining sites, with the exception of Site 1MD4, which showed extremely high biomass at the 0 m station, biomass underneath the cages was similar to the levels reported from the reference stations and peaked, depending on the site, at the 10 or 25 m station.In terms of community structure, (Table 5), the 6 sites showed a pattern of spatial change with distance from the cages.At most of the sites (1MD1, 1MD4, 1MD5 and 1MD6) at stations up to 10 m distance from the cages, the macrofauna was dominated by Capitella cf.capitata, which comprised, on average, more than 44% of the total abundance.In particular at Site 1MD6, stations up to 10 m distance were almost exclusively dominated by Capitella cf.capitata, which comprised more than 85% of the total abundance.species, namely Malacoceros fulginosus and Neanthes caudata, also contributed substantially to the total abundance up to 10 m distance from the cages at most of the sites (15 to 46%).At Sites 1MD2 and 1MD3, Capitella cf.capitata was present in proximity to the cages but in low proportions, whereas Protodorvillea kefersteini comprised more than 16% of the total abundance.Finally, at stations close to the cages in Site 1MD2, Oligochaeta contributed significantly to the total abundance representing, on average, more than 34%.At the other end of the enrichment gradient, the polychaete species Monticellina heterochaeta, Levinsenia gracilis, Cirrophorus branchiatus, Cossura costata and Magelona alleni were consistently found at the reference stations and at distances greater than 25 m, but never near or under the cages.The Shannon-Wiener diversity index and species richness (Fig. 5c,d) increased at all sites with increasing distance from the cages and obtained maximum values at the reference stations.The only excep tion to the above pattern was Site 1MD6, for which the maximum number of species was found at the 50 m station.

Flux and benthic impact indicators
The observations for the sites fell into 2 categories: Category 1 sites where macrofaunal abundance increased to a peak near the cages, followed by a decrease with increasing distance (1MD1, 1MD2, 1MD3, 1MD6), and Category 2 sites where a peak did not occur (1MD4, 1MD5).For Category 1 sites, reference stations had more than 96 species and generalised relationships were formed for these sites.Category 2 sites did not fit the generalised relationship and are plotted separately owing to the differing trends in abundance and low species numbers.Relationships were established between predicted flux and species, abundance and biomass (Fig. 6), and abundance:species ratio (A/S index), Shannon-Wiener index (H ', log 2 ), Simpson's index (1 − lambda) and BFI (Fig. 7).No significant relationships were established between predicted flux and evenness (J '), percent carbon, percent nitrogen, chlorophyll a, and phaeopigments.As expected, Margalef's species richness index showed a similar relationship to species, so is not presented.There is a gap in observational data for model predictions between 1272 and 3140 g m −2 yr −1 between the 10 and 25 m stations.
The model's capability for predicting species number at Category 1 sites (Fig. 7) was classed as 'Excellent' (Table 4) (R 2 = 0.81), with species number increasing from fewer than 20 (flux ≈ 10 000 g m −2 yr −1 ) to more than 90 for less impacted stations (flux < 2000 g m −2 yr −1 ).For Category 2 sites, 1MD4 had low numbers of species (55) at the reference station, and 1MD5 had more than 160 species at the reference station.Generalised abundance curves were classed for each mean are shown as 'Excellent' (R 2 = 0.56) (Fig. 6).Abundance peaks at 3 sites, 1MD1, 1MD2 and 1MD6, were at approximately the same modelled flux (4545 to 5501 g m −2 yr −1 ), whereas the fourth site (1MD3) showed a peak in abundance at a station with lower predicted flux (1272 g m −2 yr −1 ).MD4 and MD5 did not exhibit abundance peaks.Peaks in biomass were found at all sites, but at different locations on the flux gradient, which contributed to considerable scatter (model classed as 'Good', R 2 = 0.16) (Fig. 6).Generally, biomass was low at high-flux stations, increased to a peak and then decreased at increasing distance from the farm.For 1MD2 and 1MD3, peaks in biomass also coincided with peaks in abundance, but these peaks occurred at different locations on the modelled flux gradient (5284 and 1272 g m −2 yr −1 , respectively).For Sites 1MD1 and 1MD6, peaks in biomass and peaks in abundance did not occur at the same station.
Three of these peaks in biomass (1MD1, 1MD3 and 1MD6) occurred at relatively low modelled flux stations (<1500 g m −2 yr −1 ), and the fourth (1MD2) occurred close to high-flux areas (5284 g m −2 yr −1 ).The A/S index was low at stations more distant from the cage and at reference stations where low flux was predicted, and generally increased to a peak in A/S at stations in moderate flux areas (model classed as 'Excellent', R 2 = 0.63) (Fig. 7).In very high-flux areas under the cages, A/S is low.Of the sites where a peak in abundance was observed, the A/S ratio peaked between 3996 and 5501 g m −2 yr −1 for 3 sites (1MD1, 1MD3, 1MD6) and at 7438 g m −2 yr −1 for the fourth site (1MD2).The generalised relationship shows a maximum was present at approximately 5500 g m −2 yr −1 coinciding with peaks in abundance observed in the moderate flux zone (i.e. the zone near the farm but not immediately under the cages, with a flux of solids greater than that at the low-flux zone outside the predicted depositional footprint).Predictive capability for the Shannon-Wiener diversity index was classed as 'Excellent' (Fig. 7, R 2 = 0.67), and it increased from the high-flux, near-cage stations, to the low-flux reference and distant stations.The generalised relationship for Simpson's index increased as flux predictions decreased towards distant and reference stations (Fig. 7).The model did not predict the exceedingly low Simpson's index values at high-flux zones for 1MD6, and this resulted in the low R 2 of 0.16.Redox decreased towards the high-flux stations under the cages, but there is a high degree of scatter in all predicted flux zones for this indicator (Fig. 7, R 2 = 0.37).

Flux and indicator species and families
Relationships were established between modelled flux and relative abundance and presence/absence of indicator species and families.Using an absence criterion of relative abundance of <1%, all of the pollutant intolerant indicator species were absent at stations receiving greater than 1500 g m −2 yr −1 TDS (100% of data fit this relationship).By using the absence criterion of relative abundance of < 5%, pollutant tolerant indicator species were determined to be absent at stations receiving less than 1500 g m −2 yr −1 (86% data fit).Pollutant intolerant families were determined to be absent above a predicted flux of 1500 g m −2 yr −1 when a criterion of absence of <1% relative abundance (100% data fit) was used.By using the same criterion, pollutant tolerant families were determined to be absent at stations with a predicted flux less than 1500 g m −2 yr −1 (75% data fit).

Model capability: flux
In the high-flux zones underneath the cages in the 24 h experiments, the model predicted flux at the seabed more accurately than flux directly underneath the net or in the water column.For longer deployment periods and up to 50 m from the cages, the model predicted low (≤ 500 AFDW g m −2 yr −1 ) and high flux (> 2500 AFDW g m −2 yr −1 ) better than mid flux (500 to 2500 AFDW g m −2 yr −1 ).For these longer deployments, subtraction of background sedimentation in low-and mid-flux zones and replication of sediment traps to determine background sedimentation proved essential.
The model could not predict flux variability with traps spaced 2 m apart in the high-flux zone.However, it is not necessary for the model to resolve differences in deposition on this spatial scale to predict longer-term benthic impacts along a transect where stations are normally at least 5 m apart.Benthos within the high-flux zone (i.e.underneath the cages and at the cage edge) is normally sampled at a single station, so it is not essential for the model to resolve differences in flux at locations 2 m apart.This result emphasises the importance of adequate replication in benthic monitoring and sediment trap studies to capture spatial heterogeneity.Differences in observed flux between closely spaced traps are likely to be caused by particles exiting the cage domain at varying locations, which is affected by the shape of the net bottom, the degree of fouling, where the farmer stands when feeding the fish and fish activity and behaviour.Examination of the NET sediment trap deployments shows that in 2 of the 3 experiments, the trap located underneath the cage centre had much higher flux than adjacent traps (by factors of 2.3 and 1.4).This indicates a higher proportion of material is exiting the cage at the centre of the net bottom, which could explain the poor model performance in the NET experiments.
Excess biofouling of some nets resulted in a source of biodeposits into the NET traps in one experiment, and we could not easily separate these from waste feed and faeces.Stable isotope analysis, which quan-  tifies δ 13 C and δ 15 N concentrations in trap contents, can verify the extent and severity of deposition footprints from fish farms (Sarà et al. 2004, Yokoyama et al. 2006, Holmer et al. 2007).It is possible that this technique could help to distinguish fish feed and faeces from biofouling organisms if these are sufficiently different in stable isotope profile.Biofouling can increase biodeposition to sediments underneath shellfish farms, resulting in an additional source of material to sediment traps (McKindsey et al. 2009, Weise et al. 2009).The present study chose experimental cages with minimal fouling, but where biofouled experimental cages cannot be avoided, attempts should be made to remove non-fish wastes from trap contents if possible, or to add a source term in the model to account for this additional flux.
The accurate determination of trap and cage positions in the model is important.In the high-flux zone, model performance was defined as poor for the WC traps, but the position of the WC traps could not be exactly defined in the model grid, as these were suspended on a line in the water column from one side of the cage to the other.Although the depth of the traps were checked by divers, the exact position of the traps in the model grid was ±1.5 m in both the horizontal and vertical planes, and this could have led to inaccuracies in model predictions.By contrast, however, the position of the SB traps was measured accurately by using a weighted rope with measured distances.
Model performance was improved when the effects of wild fish feeding on pellets were included in the simulations.The improvement in model performance when the effect of wild fish is included is shown in Fig. 2, and is consistent with observations of fish behaviour at the site.Although the present study did not quantify the percentage of uneaten feed or the percentage of uneaten feed that was removed by wild fish in the field for each site, optimisation of model performance resulted in para meters of feed wastage of 5%, and 50% of this was removed by wild fish feeding.These factors will vary between farms as a result of different husbandry practices, and the Modelled flux + 1 (g m -2 yr -1 ) 5 1 10 1 10 2 10 3 10 4 10 5 1 10 1 10 2 10 3 10 4 10 5 1 10 1 10 2 10 3 10 4 10 5 1 10 1 10 2 10 3 10 4 10 5 Fig. 7. Benthic response model for Category 1 generalised relationships (1MD1, 1MD2, 1MD3, 1MD6) and Category 2 data (1MD4, 1MD5) for A/S, Shannon-Wiener index (H ', log 2 ), Simpson's index and biomass fractionation index (BFI).Category 1 sites showed a peak in macrofaunal abundance whereas Category 2 sites did not influence of wild fish feeding will vary between seasons.The contribution of a high wastage rate of feed to seabed flux in the summer months when feed input is high will be reduced by a higher level of wild fish feeding in these months.Conversely, during winter months, the effect of wild fish feeding on uneaten pellets is likely to be lower than in the summer months when these validation experiments were undertaken.The influence of uneaten pellets on the predictions of sediment trap flux and benthic response decreases with increasing distance from the farm, as the proportion of feed to faecal particles decreases.A feed particle settling at 10.8 cm s −1 in 25 m depth will deposit on the bed at 5 m from the cages (assuming a mean current speed of 2 cm s −1 ).Thus, these parameters relating to feed and wild fish primarily influence the flux to the traps directly underneath the cages and benthic sampling stations at 5 m from the cages.At more dispersive and deeper sites than those studied here, waste feed parameters will influence predictions at more distant stations.
Several studies (e.g.see 'Introduction') have shown the importance of marine fish farms as attractors of wild fish, and based on a variety of techniques (e.g.lipid analysis, Fernandez-Jover et al. 2007) to identify wild fish, at least one of the reasons for this attraction is to feed on fish wastes.However, estimating the scale of the trophic transfer of fish wastes to wild fish is very difficult experimentally as wild fish exist in a relatively open system.Our sediment trap survey at 3MD8 was, to our knowledge, the most extensive ever carried out at a fish farm, and fitting our particle tracking model to optimise the relationship between predictions and observations of waste particle flux probably gives a reasonable indication of the scale of this trophic transfer.If this result is typical of fish farms in the oligotrophic eastern Mediterranean Sea, then it may help to better understand the relationship between fish farming and the structure and dynamics of this ecosystem.

Model capability: benthic response
A benthic response model was validated adequately, and generalised relationships were established for species, abundance, A/S, Shannon-Wiener H ' and Simpson's index.Less useful relationships were established for biomass and redox.The relationships were established for sites where an abundance peak was observed at near-cage stations (Category 1 sites).Two other sites did not fit these generalised relationships (Category 2 sites), but nevertheless showed expected trends for species, A/S and Shannon-Wiener and Simpson's indices.
An absence of flux predictions between 1272 and 3140 g m −2 yr −1 was caused by predictions at 10 m stations being above 3140 g m −2 yr −1 and predictions at 25 m stations being below 1272 g m −2 yr −1 .This interesting result suggests that consistently sampling at 10 and 25 m stations at all sites resulted in undersampling of a benthic community exposed fluxes between 1272 and 3140 g m −2 yr −1 .This similarity occurred across all sites, even though all sites were modelled with different environmental and husbandry characteristics.However, benthic response was not tested at these fluxes, but the accuracy of flux predictions was tested in the trap studies, where model performance in the mid-flux zone was classed as 'Good' for 27 traps.
The maximum distance that benthic impact was observed in this study was 25 m from the farm, which was supported by observations of low-flux in sediment traps at this distance.These results do not support the findings of some other studies, where effects of the farm are found up to a distance of 1000 m (Sarà et al. 2004).Theoretical consideration of particle trajectories further suggests that for the measurements of current speed, depth and particle settling velocity, significant deposition at a distance of 1000 m is highly unlikely.Such deposition is likely to occur from wild fish faeces or by redeposition from resuspension.The latter, however, is unlikely given that the critical erosion threshold for fish farm wastes is typically around 10 cm s −1 and a very small percentage of bed current speed exceeded this threshold for these study sites.As offshore aquaculture is developed in deeper and more dispersive sites, the areal extent of deposition footprints will increase while the intensity of effect will decrease.
Increasing TOC levels are correlated with decreasing species diversity at fish farms (Kalantzi & Karakassis 2006), and several studies have considered the relationship between macro benthic community response and sedimentary organic carbon.For example, Hyland et al. (2005) suggested that risks of species richness reduction are low at TOC < 1%, high at TOC > 3.5% and intermediate between these thresholds.However, predictions of TOC from TDS flux are difficult to make as a large number of factors may intervene between TDS deposition and carbon accumulation in sediments.The present results show that a modelled flux rate of 1500 g m −2 yr −1 TDS is a watershed between those in which the pollution tolerant and intolerant taxa chosen would typically be found.Cromey et al. (2002) presented data for sam-ples taken from 2 Scottish fish farms on the relative abundance of Trophic Groups 1 to 4, as listed in the Infaunal Trophic Index (ITI, Word 1980); these groups are considered to have increasing tolerance to pollution in that order.Plotting these salmon farm data against modelled prediction on TDS flux showed that taxa in Groups 1 to 3 were present only in low relative abundance at fluxes above 1000 g m −2 yr −1 .Group 4 taxa were dominant (> 80% relative abundance) at fluxes above 1000 g m −2 yr −1 and were also present (0 to 80% relative abundance) at stations with modelled fluxes of 0 to 1000 g m −2 yr −1 .The nomogram of Hargrave et al. (2008) shows that at predicted TDS fluxes between 730 and 1825 g m −2 yr −1 , indicator species diversity changes from moderate to reduced.Although not directly comparable, the results presented here are consistent with other evidence that describes the transition between polluted and unpolluted sediments.

Model limitations
The model should be used to predict flux and benthic response with special regard to the model accuracy specified.Model accuracy varies depending on the magnitude of flux and each benthic impact indicator.Using detailed husbandry data in the model is important, and the use of monthly summarised husbandry data for the whole farm will be limiting owing to the range of fish size and species being farmed within a cage group.The model should not be used to predict differences in flux or benthic impact for stations less than 5 m apart; averaging of observed flux underneath cages is recommended before comparing with predictions in validation exercises.This model has not been tested in areas of hard substrate where there are underwater cliffs, nor has it been tested where there are significant resuspension effects in Mediterranean waters.The model does not include flocculation or disaggregation behaviour of particles.

CONCLUSIONS
We have developed a model, MERAMOD, of the benthic impact of sea bass and sea bream aquaculture in the eastern Mediterranean Sea and calibrated it with a comprehensive series of sediment trap deployments that include multiple short deployments at fine resolution in the near-field around fish farms and a major 57 trap study in the 5 to 50 m radius around one fish farm over 13 d.Two of the farms studied did not have peaks in macrofaunal biomass and abundance but had similar relationships between predicted flux and various benthic indices as the other 4 farms studied.The model includes the facility to remove a percentage of the waste feed that is eaten by wild fish and therefore does not contribute significantly to local benthic impact.MERAMOD is also highly adaptable for farms with co-culture of species with faecal material that has different settling velocities and where several size classes of fish (cohorts) are present on the farm.The model provides a sound basis for regulating, planning and monitoring the large fish farming industry in the eastern Mediterranean and has the capability to be adapted to other environmental and husbandry settings.

Fig. 2 .
Fig. 2. Observed and modelled ash-free dry weight (AFDW, n = 57) for Site 3MD8 showing the improvement in model performance by scenario.(a) Scenario a includes background sedimentation and (b-d) Scenarios b, c and d have background sedimentation subtracted from observations in which 0, 100 and 50% (optimised) of waste feed pellets were removed by wild fish, respectively

167Fig. 3 .
Fig. 3.The 24 h sediment trap experiments showing mean (± SD) observed and predicted total dry solids (TDS) fluxes (total n = 72 traps, n = 12 traps per experiment): (a) NET traps deployed attached to nets, (b) SB traps on seabed and (c) WC traps suspended in the water column.The factor by which the model is different to the observed TDS and the number of samples (n) for each mean are shown

Fig. 4 .
Fig. 4. Geochemical parameters measured in the sediment at the reference stations (Ref) and at various distances (m) from the cages.(a) Median diameter (MD) of the sediment; (b) redox potential (Eh) at 4 cm depth; (c) total organic carbon (TOC)

Table 2 .
Summary of hydrographic and wind data at 7 sites showing deployment time (D) of sediment traps in long-term studies on benthic response studies and short-term particle track studies.Sites with prefix 1 were measured in July 2001 except 1MD2(December 1998) and 1MD4 (July 1996).Sites with prefixes with 2 and 3 were measured in March and October 2002, respectively.A and B: experiments conducted at the same site and in the same month on different days.Meter locations: S, surface; M, mid-water (i.e.net bottom); B, near-bed.Meter type: RCM7 = rotary; S4 = electromagnetic; ADCP = Acoustic Doppler Current Profiler.Current speed was measured as cm s -1 , wind speed (in italics) as m s -1 .See 'Bathymetry and hydrography' for sampling intervals of wind speed measurements.°T: degrees true.Dispersion coefficients were measured for deployment lengths of 1 to 7 h

Table 3 .
Two other polychaete Summary of model accuracy and the range of predicted flux values for which the model was tested for during all sediment trap experiments.The MD8 experiment was 13 d in length.The mean difference between observed and predicted flux for 3MD8 (mid-flux zone) was ± 29% of the predicted value.The 24 h statistics are combined from studies

Table 4 .
Statistical tests using the linear regression equation Y i = β 0 + β 1 X i for observed (Y i ) and predicted (X i ) values for the benthic response model and sediment trap studies.Scenario a for Site 3MD8 includes background sedimentation and Scenarios b, c and d have background sedimentation subtracted from observations in which 0, 100 and 50% (optimised) of waste feed pellets were removed by wild fish, respectively.Scenario d-Low, -Mid, -High apply to low-, mit-and high-flux zones, respectively.BFI: biomass fractionation index; NET: attached to the net bottom; SB: on the seabed; WC: suspended in the water column.See 'Materials and methods: Mass balance checks ...' for description of class designations

Table 5 .
Average relative abundance (%) of macrofaunal species comprising more than 5% of the total abundance of macrofauna found at least at one station over all sites (+ indicates presence of <1%; blank indicates absence).Species are arranged by decreasing average abundance at the stations in proximity to the cages (at 0, 5 and