- Email: [email protected]
Impact of aquacultures on the marine ecosystem: Modelling benthic carbon loading over variable depth
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/ecolmodel
Impact of aquacultures on the marine ecosystem: Modelling benthic carbon loading over variable depth Marko Jusup ∗ , Sunˇcana Geˇcek, Tarzan Legovi´c – Boˇskovi´c Institute, Bijeniˇcka Cesta 54, 10002 Zagreb, Croatia Ruder
a r t i c l e
i n f o
a b s t r a c t
Article history:
In order to be used within Environmental Impact Assessment study, we have developed a
Received 2 August 2005
three-dimensional particle tracking model for prediction of benthic carbon loading (BCL)
Received in revised form
caused by fish farms. The model is based on stochastic differential equations for particle
15 August 2006
transport consistent with the well-known semi-empirical advection/diffusion equation. It
Accepted 24 August 2006
requires only easily obtainable input data in the form of measured current record, the source
Published on line 5 October 2006
location and a specification of local bathymetry. The model accounts for advection by longterm residual and tidal currents, turbulent diffusion, realistic bathymetry and variations in
Keywords: Fish farms
daily (monthly or yearly) emissions from fish farm. Here, we concentrate on the changes in sedimentation pattern caused by various bathy-
Benthic carbon loading
metric shapes. Examination of idealized cases reveals where and why we can expect the
Model
worst impact on benthic communities. For future reference, these results will be included
Variable depth
into guidelines for fish farming. © 2006 Elsevier B.V. All rights reserved.
1.
Introduction
Mariculture has already become very important area of food production all over the world. Potentially high profits have attracted investors in the last two decades triggering its rapid development. At the same time, many concerns have been raised about possible negative environmental impact, especially at locations where regional economy is tourism based. Main sources of environmental problems in connection with fish farming are enrichment of both water column and benthos with organic material, and physical and biological disturbances caused by cage structures (Silvert and Sowles, 1996). Therefore, in order to ensure the long-term use of farming sites Ervik et al. (1997) proposed three basic sustainability requirements and later Hansen et al. (2001) presented details of appropriate monitoring program. In these papers most concern is devoted to prevention of disappearance of benthic
∗
infauna due to excessive accumulation of particulate organic wastes. The first step toward resolution of this issue is development of a model that accurately estimates benthic nutrient loading. Hence, we present a tool, named KK3D, for prediction of dispersal of particulate organic matter from fish farms. It belongs to a broad class of Lagrangian stochastic models (LSM) which have been widely used in atmospheric sciences (see Wilson and Sawford, 1996, for a review). In oceanography the interest for them is relatively low and some of the probable reasons are discussed in Section 3. The model is based on stochastic differential equations consistent with the well-known semi-empirical advection/diffusion equation. It requires only easily obtainable input data in the form of measured current record, location of the source and the specification of local bathymetry. The model accounts for advection by long-term residual and tidal currents, turbulent diffusion, realistic bathymetry and variations in daily (monthly
Corresponding author. Tel.: +385 1 4680230. ˇ ´ E-mail addresses: [email protected] (M. Jusup), [email protected] (S. Gecek), [email protected] (T. Legovic). 0304-3800/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2006.08.007
460
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
or yearly) emissions from fish farms. Special attention is given to changes in the sedimentation pattern due to variable bathymetry. The paper is organized as follows. In Section 2 an overview of other benthic nutrient loading models is given. It is followed by a brief review of LSM theory and details of carried simulations in Section 3. Section 4 displays the results. The paper is concluded with a discussion of results in Section 5.
2.
Overview of other models
Modelling offers methods for analyzing various aspects of aquaculture like: - identification of suitable sites by GIS technique (Congleton et al., 1999), - assessment of carrying capacity for multispecies culture by coupling physical and biogeochemical models (Duarte et al., 2003), - simulation of the dynamics of an integrated experimental aquaculture plant (Buonomo et al., 2005), long before the production is initiated. Environmental impact estimation of marine fish farms, on both local and regional level, is no exception. So far, much attention was given to particle tracking simulations that predict distribution of particulate organic wastes on the seabed (i.e. the sedimentation pattern or the deposition footprint). One of the simpler approaches (Gowen et al., 1994) is based on the observation that the dispersion distance D is related to current speed U, water depth H and settling velocity of pellets w by D=
UH . w
(1)
It represents the distance that particles travel from the cage before settling on the bottom. The presence of water depth in Eq. (1) makes the “Gowen” method suitable to include variable bathymetry. However, it was Hevia et al. (1996) who first extended the model to account for bottom topography and vertical variations in current speed. The importance of turbulent diffusion in horizontal dispersion was also realized, resulting in an additional term in the relation for dispersion distance (Gillibrand et al., 2002) UH D= + w
2KH . w
where K is the turbulent diffusion coefficient. Stigebrandt et al. (2004) adopted somewhat different approach. They estimate the current speed variance 2 from the measurements and use it to determine the dispersion length T = H/w of a given site. For each dispersion length T there is a pre-calculated, dimensionless, normalized loading (sedimentation) function = (r) which relates emission from the net pens F1 with benthic carbon loading F2 through the relationship F2 (r) = (r)F1 . The distance from the cage center is denoted by r. Note how this procedure should be computationally inexpensive once
the loading functions have been determined. However, that would be an impossible task without the fact that current velocities appear to be approximately normally distributed (Green and Stigebrandt, 2003). DEPOMOD, developed by Cromey et al. (2002a), is a tool based on direct particle tracking. It was validated in sediment trap experiments and showed good agreement with field data with accuracy of ±20 and ±13% for a dispersive and depositional site, respectively. The model successfully links solids accumulation to benthic indices and total individual abundance. Two processes determine the path of a particle: advection by measured currents and turbulent diffusion modelled with random walk. This approach is analogous to the methods used in KK3D. DEPOMOD is accompanied by the validated resuspension model (Cromey et al., 2002b) which is a rare occurrence within the field. All the mentioned particle tracking models are characterized by physical simplicity and reasonable predictive capability. However, none of them is equivalent to the advection/diffusion equation which governs the tracer dispersion from a given source. The use of LSM enabled us to achieve that goal. Consequently, KK3D is a tool based on proven physical principles and underlying stochastic differential equations are quite simple to implement as a numerical simulation on a PC. Furthermore, recent inclusion of the adsorption term naturally accounted for nutrient loss caused by leaching. Information on the chemical composition of organic matter in faeces reaching the bottom allows the calculation of oxygen demand (g O2 /g C). For that purpose stoichiometric analysis similar to Krasakopoulou et al. (2005) was adopted. Thus, KK3D addresses the chemical and biological processes which determine the final fate of nutrients emitted from a fish farm. However, full description of these additions is beyond the scope of this paper. We cannot claim that KK3D resolves all questions regarding the sedimentation pattern. The model user should be aware of an important limitation because there is no resuspension of sediment. It is a physical process which, beside the mentioned chemical and biological processes, affects the final distribution of organic waste on the seabed. Resuspension is usually split into four compartments: erosion, transport, deposition and consolidation. Erosion occurs when the shear stress on the bottom = U∗2 becomes greater than the critical erosion shear stress ce . Here, is the density of seawater and U* is the bed shear velocity which is estimated from the measured current speed near the bottom Uz assuming a logarithmic profile U∗ =
Uz . ln(z/z0 )
Thus, the hydraulic roughness length z0 , as an important characteristic of the sediment type, plays a role in the formation of sedimentation pattern. Transport is modelled by the same equations used for determining the initial distribution of waste material. Deposition occurs when the shear stress on the bottom becomes lower than the critical deposition shear stress cd . It is assumed that cd < ce so that erosion and deposition are mutually exclusive. Consolidation is a function of
461
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
time spent on the seabed. The longer the particle stays on the bottom, probability of being eroded approaches zero. To conclude, only a full resuspension model can account for properties of particular sediment type and this has not yet been included into the KK3D tool. Finally, we fell that the oxygen demand should be related to the flux of oxygen into the sediment and oxygen consumption by the microbial degradation. A promising attempt to resolve that problem was made by Stigebrandt et al. (2004).
3.
Methods
It is a common practice to breakdown the sea current velocity into two parts: (1) slowly varying component due to tidal, thermohaline and wind forcing and (2) random component related to turbulence that quickly changes in time. Former part is responsible for displacement of waste material from the source (advection) while the latter part causes the dispersion (turbulent diffusion). Therefore, it is tempting to adopt the semi-empirical approach (Monin and Yaglom, 1971) and model the transport of any passive tracer emitted from an aquaculture by using the advection/diffusion equation: ∂C ∂C ∂ = + vi (x, t) ∂t ∂xi ∂xi
Kij (x, t)
∂C ∂xj
,
vi (x, t) + vz ıi3 +
∂Kij (x, t) ∂xj
dt +
vi C = Kij
∂C . ∂vj
Particle tracking model described above requires two types of input data:
(2)
where C is a tracer concentration, vi , i = 1, 2, 3 the current velocities and Kij , i, j = 1, 2, 3 are the eddy diffusivity coefficients. This is an Eulerian approach and we can replace it by a Lagrangian technique called Random Displacement Model (RDM) which is consistent with Eq. (2) (Rodean, 1996). However, in order to predict the distribution of particulate organic material on the seabed we must calculate the paths of heavy particles and not the passive tracer. According to Wilson (2000), RDM easily resolves this problem by superimposing to it the particle settling velocity vz . Thus, the position of a particle in RDM is given by stochastic differential equations:
dxi =
2005). Furthermore, there is a problem of so-called trajectory crossing. Initially the heavy particle position will coincide with the position of some fluid element. As the time passes, particle drifts away from the fluid element because of its settling rate, leaving the eddy in which correlation with initial velocity still exists. Therefore, it is suggested to reduce the heavy particle correlation timescale with respect to the fluid element timescale (e.g. Wilson, 2000), but this is a heuristic approach. There is also a possibility that due to finite inertia, heavy particle is unable to follow rapid velocity fluctuations, resulting in a different velocity variance compared to the fluid element. Our objective was to construct a model that is usable in many different locations prior to installation of an aquaculture. In such situations, at best, measured current profiles are available as input data. So to avoid unnecessary complexities and account for dispersion, one must impose reasonable assumptions about joint moments of concentration and velocity fluctuations. The simplest and widely used example is the eddy diffusion hypothesis
2Kij (x, t)dWj ,
(3)
where dWi , i = 1, 2, 3 are the standard Wiener processes. Sedimentation pattern is expressed in form of benthic carbon loading (BCL) in gC/m2 day. Wilson (2000) classifies RDM as a “zeroth-order” LSM noting its main weakness; if xi (t) is treated as Markovian random variable, the model cannot be valid close to the source, that is for travel times shorter than the velocity correlation timescale. The more advanced, “first-order” LSM, called Langevin Equation Models (LEM), have been used extensively in atmospheric research owing much of their popularity to analytical description of turbulence properties in the surface boundary layer. For example, Kurbanmuradov and Sabelfeld (2000) have shown how a number of simplifying assumptions lead to unique selection of three-dimensional LEM for turbulent diffusion in the atmosphere. This turns out to be an impossible task in oceanographic application (Brickman and Smith, 2002) where turbulent quantities are provided by circulation models with incorporated appropriate closure scheme (e.g. Warner et al.,
1. Measured current record from a representative location must be obtained. Slowly varying tidal components as well as long-term residual mean should be extracted from such a record while short-period fluctuations related to turbulence are of no interest here since we adopted the eddy diffusion assumption. Processing of current meter data by classical harmonic analysis (Pawlowicz et al., 2002) provides satisfactory results returning the following information about each tidal component: - is the frequency of tidal component in cycles per hour (cph), - a is the major semi-axis of tidal current ellipse (cm/s), - b is the minor semi-axis of tidal current ellipse (cm/s), - is the inclination of tidal current ellipse, - g is the Greenwich phase, - V(t0 ) is the astronomical argument at the moment t0 , - u is the nodal phase modulation, - f is the nodal amplitude correction. From these parameters, the model generates tidal current according to following relations (Foreman, 1978): a+ =
a+b f, 2
a− =
a−b f 2
ε+ = − g + (u + V(t0 )),
ε− = + g − (u + V(t0 ))
(4a) (4b)
v1 (t) = a+ cos(ε+ + 2 t) + a− cos(ε− − 2 t)
(4c)
v2 (t) = a+ sin(ε+ + 2 t) + a− sin(ε− − 2 t).
(4d)
Thanks to Eqs. (4a)–(4d) KK3D is capable of truly forecasting the sedimentation pattern and not only hind-casting as it is the case when measured current record is directly loaded into
462
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
the model. Extreme events, like storms, can be modelled by modifying the residual current for a short period of time. Modification should be made according to measurements or results from an appropriate hydrodynamic model. Also, activities like net cleaning are easily included into the prediction. The user should simply increase farm output during the desired time period. 2. Rectangular grid containing bottom topography is also needed. Grid cell resolution is arbitrary and depends only on available data from a location of interest. Cromey et al. (2002a) used the resolution of 10 or 25 m depending on the size of the sedimentation pattern and in realistic application one can hardly expect higher precision due to instrument limitations. Here, we have concentrated on emissions from a single net pen, which cover relatively small area on the seabed. Therefore, slightly denser grid is used for obtaining the results discussed in this paper. To start particle tracking we also need the magnitude of settling velocity vz . The velocity is chosen randomly for each particle from the normal distribution with the mean = 3.2 cm/s and the standard deviation = 1.1 cm/s (Cromey et al., 2002a). Approximate solutions to Eq. (3) are obtained by means of the Euler–Maruyama method or, as Kloeden and Platen (1992) put it, an Ito–Taylor scheme of strong order 0.5:
xi (t + t)xi (t) +
vi (x, t) + vz ıi3 +
∂Kij (x, t) ∂xj
t +
2Kij (x, t) t j .
Here, t is the time-step and j ∈ N(0, 1). In the atmospheric surface boundary layer, it is convenient to choose coordinate axes so that the x1 -axis coincides with the direction of the mean wind velocity. Then, the crosswind velocity fluctuations are symmetrical with respect to the planes x1 = 0 and x3 = 0, and we have v1 v2 = 0 and v2 v3 = 0. In addition, since the Ox1 , Ox2 and Ox3 are clearly the preferred axes of the flow they are often considered as principal axes of the diffusivity tensor. Within the sea, the direction of mean flow can change with depth making the above reasoning inapplicable. However, Monin and Yaglom (1971) note that “all qualitative predictions of semi-empirical diffusion theory can hardly depend on the presence or absence of terms with coefficients K13 and K31 ” (authors considered the inclusion of these terms into Eq. (2)). Accordingly, our model includes only diagonal elements of the diffusivity tensor. Besides, these coefficients are depth dependent to account for the existence of surface and bottom boundary layers. Their mean values correspond to the values observed in Scottish inshore waters (K11 = K22 = 0.1 m2 /s and K33 = 0.001 m2 /s) and are used for modeling by Gillibrand et al. (2002) and Cromey et al. (2002a). Thus, based on published results, we expect that the key properties of sedimentation pattern such as its shape and size, the location of maximum BCL, etc., can be estimated with acceptable accuracy. Initial particle position in x1 Ox2 plane is chosen randomly from a uniform distribution assuming that the fish can be located anywhere within the cage limits immediately prior to defecation. In all simulations the cage radius was set to 20 m and the center was placed at x1 = x2 = 0. Let us note that fix-
Table 1 – Tidal ellipse parameters for constituents used in all simulations Tidal constituent
O1 K1 M2 S2
Frequency (cph) 0.0387307 0.0417807 0.0805114 0.0833333
Major semi-axis (cm/s) 2.00 5.32 8.50 5.20
Minor semi-axis (cm/s) −0.23 −0.36 −0.32 −0.27
Negative minor semi-axis indicates clockwise rotation of tidal current phasor. Inclination was set to zero. Note the highly polarized tidal currents typical for the Adriatic Sea.
ing the starting position (for example, to the cage center) does not affect the sedimentation pattern at all in sufficiently deep water. In such conditions, the deposition area is much larger than the cage area and the exact initial position of the particle is hardly relevant. However, in shallow water fixing the starting position can result with under-prediction of sedimentation pattern size. The model is far more sensitive to the choice of initial vertical coordinate x3 because dispersion is depth dependent, i.e. deposition area increases rather quickly with depth. Since there is observational evidence that defecation occurs near the surface immediately after feeding, we have chosen x3 = −5 m for all particles. The number of tracked particles in each simulation was set to 2 × 105 since no changes in results occurred with higher values. Output data are passed to visualization module, which produces contour plots revealing details of predicted sedimentation pattern. Unless we simulate exceedingly large number of sample paths, the stochastic nature of employed algorithm will cause appearance of small patches in final images. Clearly, these patches are insignificant as they disappear with large number of tracked particles. To avoid spending too much computational time we recommend the removal of model patchiness by a two-dimensional wavelet analysis. In this study, dispersion is simulated over 1-day period, but that value can easily be changed in order to meet specific requirements. For example, sometimes it is useful to predict dispersion when strong currents towards the shore dominate the aquaculture area. The model also accounts for uneven emission rate of waste material during the simulated period. This is necessary as we expect heavier impact in direction of prevailing current when emission rate reaches maximum. Here, we have concentrated on changes in sedimentation pattern induced by various bathymetric features. Consequently, for all examined cases the residual and tidal currents are the same. In Table 1 information on used tidal constituents is given. Values are typical for the Adriatic Sea. One should keep in mind that all studied cases are idealized.
4.
Results
In order to investigate the validity of described model we have run a series of simulations that were supposed to satisfy the following criteria: (i) in case of disabled advection, the model should reproduce the results of Stigebrandt et al. (2004) illustrated in dispersion submodel section of their paper and (ii) in
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
463
Fig. 1 – Schematic view of bathymetry over a transect y = 0 for all simulations: (a) flat seabed, (b) 2.5◦ slope, (c) 5◦ slope, (d) 5 m deep pit centered at x = −100 m and (e) 5 m high hill centered at x = −60 m. Arrows in (a) indicate direction and magnitude of residual current.
case of bottom slope, the model should reproduce the effects reported by Hevia et al. (1996). Both criteria were successfully fulfilled; instance (i) is not further discussed here because it is not directly related to the main topic, but (ii) is explained in detail below. Computer simulation results are illustrated in Figs. 1 and 2. Fig. 1 shows the crosssectional view of bathymetry along the line y = 0 as well as the position of net pen for all studied cases. Contour plots in Fig. 2 display the corresponding sedimentation patterns. The first simulation was made for the flat seabed and it is considered a reference point to which all other cases are compared. Before examining the changes in sedimentation pattern caused by several different bottom topographies, we must introduce few important concepts. Let us denote the maximum BCL by m so that the values of BCL fall within the segment S = [0, m]. Subdivision of S into six subintervals S1 = [0, x1 ], S2 = [x1 , x2 ], . . ., S6 = [x5 , m] (where S1 , . . ., S5 are equally long) provides the basis for contour plots shown in Fig. 2, that is, isolines (lines of constant BCL) are drawn at levels xi (i = 1, . . ., 5). Note how such procedure exposes
the details of each sedimentation pattern and allows one to define the maximum deposition region as a set of coordinates M = {(x, y) : BCL(x, y) ∈ S6 }. If surface z = z(x, y) represents the seabed (i.e. z is the depth at (x, y)), it is possible to approximate an area A corresponding to the region M. Parameter A, named the maximum deposition area, provides valuable information for comparing various test cases. Other features revealed by contour plots are the shape of M and the position where maximum BCL occurs, i.e. the coordinates (x0 , y0 ) such that m = BCL(x0 , y0 ). The impact on the bottom is considered negligible if the BCL level is below 0.1 gC/m2 day. Table 2 is a convenient summary for the presentation of results that follow below. For the reference simulation, constant depth level was set to 35 m (Fig. 1a). Sedimentation pattern is elongated in agreement with the polarized tidal current and displaced in response to residual current (Fig. 2a). Maximum deposition region is split into two parts what is, again, a consequence of polarized tidal current which reverses its direction throughout a day. With assumed emission of 10 kgC/day in the form of
464
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
Fig. 2 – BCL (in gC/m2 day) for all simulated test cases: (a) flat seabed, (b) 2.5◦ slope, (c) 5◦ slope, (d) 5 m deep pit centered at x = −100 m and (e) 5 m high hill centered at x = −60 m. Maximum BCL location is marked with ‘x’. Dashed circle is a projection of cage on the bottom.
feces, the maximum BCL reaches the value of 1.17 gC/m2 day at the distance of 84 m from the cage center. The total area affected by organic waste is 21,840 m2 and the area of maximum deposition equals 1704 m2 .
Even a small variation in bathymetry can cause major changes in the sedimentation pattern. This is shown in the simulation for 2.5◦ bottom slope (Fig. 1b). Note that this is the test case from the criterion (ii) for the validity of our model
Table 2 – A summary of the simulation results Maximum BCL (gC/m2 day) 1.17 1.55 1.92 1.20 1.56
Maximum BCL position (m) −84, 4 −72, 0 −66, −4 −114, 0 −54, 0
Total affected area (m2 )
Maximum deposition area (m2 )
21840 20828 20082 22165 21925
1704 2474 2722 1325 1810
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
as mentioned above. The depth range caused by a slope is 22–48 m with the mean of 35 m. Sedimentation pattern in Fig. 2b discloses the restricted deposition towards the shallow water because the tilt blocks the particles from reaching the coastline. The value of maximum BCL increases to 1.55 gC/m2 day and is attained at 72 m from the cage center. Total area influenced by organic material decreases to 20,828 m2 , but the maximum deposition area grows up to 2474 m2 . Doubling the bottom slope to 5◦ makes the changes in sedimentation pattern more prominent (Fig. 1c). Again, the mean depth is 35 m, but steeper slope causes wider depth range over the simulation area (9–61 m). At 150 m from the cage center towards the shallow area there is no noticeable deposition (Fig. 2c). For the reference simulation and in case of gentler slope the same happens at 210 and 174 m from the cage center, respectively. Hence, the deposition in shallow water is more restricted by higher slope steepness. The maximum BCL reaches the value 1.92 gC/m2 day and is closer to the net pen than previously (x0 = −66 m). The total area disturbed by organic waste reduces to 20,082 m2 and the maximum deposition area is 2722 m2 . Following bottom topographies were designed to illustrate how by neglecting relatively small bathymetric features one is lead to miss-prediction of all important parameters. More precisely, the model can produce wrong shape of the maximum deposition region (M) and give incorrect values for its size (A), the maximum BCL (m) and the corresponding coordinates (x0 , y0 ). Consider the 5 m deep pit with inner wall steepness of less than 10◦ (Fig. 1d). It is located on an otherwise flat bottom that is 35 m deep (as in the reference simulation). The pit is centered at 100 m from the cage center, meaning that it does not coincide with original maximum deposition region. However, the presence of such pit has important consequences on sedimentation pattern and most notable is the displacement (of one part) of M further away from the net pen (Fig. 2d). The maximum BCL of almost 1.20 gC/m2 day is now located at x0 = −114 m. At the location where maximum BCL occurs originally (x0 = −84 m), the predicted value has changed to 0.68 gC/m2 day; not even close to the highest loading. Seabed area that is worst affected is somewhat reduced (A = 1325 m2 ) as a consequence of better dispersal over the pit. Exactly the opposite changes are induced by 5 m high hill centered at the distance of 60 m from the net pen (Fig. 1e). The hill, as a projection from the bottom, prevents particles from reaching their original settling location. For this reason (a part of) the maximum deposition region is now displaced closer to the cage (Fig. 2e). Also, being the highest point on the bottom, the hill makes the dispersion less effective. As a result, the maximum BCL increases to 1.56 gC/m2 day and that value is recorded at x0 = −54 m; near the hill top. At the location where the maximum BCL occurs originally (x0 = −84 m), the predicted value has changed to 0.74 gC/m2 day. Thus, ignoring the hill we would completely miss both parameters: the highest BCL m and corresponding coordinates (x0 , y0 ). Finally, there are no major changes in the total affected area (21,925 m2 ) and maximum deposition area (1810 m2 ).
5.
465
Discussion
We have presented the capabilities of a model for prediction of BCL caused by fish farms. The model is based on strong, yet relatively simple, physical grounds and can account for number of features important in shaping the organic waste sedimentation pattern. Among these features, special attention is given to modelling the dispersion over variable bathymetry. For this purpose we have compared the results from several simulations with the sedimentation pattern obtained for the flat seabed, holding other environmental variables constant. The inevitable conclusion is that approximating the realistic bathymetry with the uniform bottom can result in serious error. It is possible to predict wrong shape of the maximum deposition region, miscalculate its area, give incorrect value and location of the maximum BCL and even fail to estimate the total area exposed to particulate waste. Thus, as others have stated before, necessity for more robust models exists, but the balance between complexity and available data from field must be retained. Accordingly, the presented model is based on data that can easily be recorded in practice and (even more importantly) with limited expenses. Modular programming concept, popular because it easily accommodates further extensions of the model (Cromey et al., 2002a; Stigebrandt et al., 2004), has been applied here in order to finally produce user-friendly application for fish farm installation and management. So far, the described BCL module has been used in several Environmental Impact Assessment studies (required by the Croatian Ministry of Environmental Protection, Physical Planning and Construction) for farms with production that exceed 50 t/year. In future we expect the model to become the standard, not only in the selection of suitable locations, but also for the improvement of farming practices. This implies achieving the maximum production that ensures long-term use of a site. For example, one way of approaching that goal is linking BCL to oxygen consumption in order to predict the change from aerobic to anaerobic conditions in the sediment. As this change may affect the farm performance, it would be useful to predetermine the optimum tonnage for the net pens (i.e. carrying capacity) at which production levels do not decrease. Future development may proceed in two directions: (1) refinement of physical aspects of the model and (2) conversion of predicted quantities to biochemical parameters. Regarding point (1) the most obvious step is constructing the module for post-depositional changes like the slumping of material down the steep slopes or resuspension of sediment into the water column. Another possibility is coupling of dispersion predictions with full hydrodynamic model that can resolve currents horizontally and calculate eddy diffusivities. However, practical use of such model is questionable if we consider the input data necessary for running it. Perhaps the most important issue of point (2) is connection of BCL with oxygen consumption in sediments as explained above. Local and regional water quality is also an important parameter. On the local level the water should have sufficient oxygen concentration and low concentrations of potentially harmful substances. Regional water quality implies insignificant contribution of mariculture to algal blooms which occur as a consequence of nutrient
466
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
enrichment. Some of these topics have already been addressed in the literature, but definitive solutions are yet to be reached.
Acknowledgments This research has been supported by the Croatian Ministry for Science, Education and Sport, the FP6 Project on ecosystem approach to sustainable aquaculture and the NorwegianCroatian Project on the particularly sensitive sea area.
references
Brickman, D., Smith, P.C., 2002. Lagrangian stochastic modelling in coastal oceanography. J. Atmos. Ocean. Tech. 19, 83–99. Buonomo, B., Falcucci, M., Hull, V., Rionero, S., 2005. A mathematical model for an integrated experimental aquaculture plant. Ecol. Model. 183, 11–28. Congleton Jr., W.R., Pearce, B.R., Parker, M.R., Beal, B.F., 1999. Mariculture siting: a GIS description of intertidal areas. Ecol. Model. 116, 63–75. Cromey, C.J., Nickell, T.D., Black, K.D., 2002a. DEPOMOD-modelling the deposition and biological effects of waste solids from marine cage farms. Aquaculture 214, 211–239. Cromey, C.J., Nickell, T.D., Black, K.D., Provost, P.G., Griffiths, C.R., 2002b. Validation of a fish farm waste resuspension model by use of a particulate tracer discharged from a point source in a coastal environment. Estuaries 25, 916–929. Duarte, P., Meneses, R., Hawkins, A.J.S., Zhu, M., Fang, J., Grant, J., 2003. Mathematical modelling to assess the carrying capacity for multi-species culture within coastal waters. Ecol. Model. 168, 109–143. Ervik, A., Hansen, P.K., Aure, J., Stigebrandt, A., Johannessen, P., Jahnsen, T., 1997. Regulating the local environmental impact of intensive marine fish farming. I. The concept of the MOM system (Modelling-Ongrowing Fish Farms-Monitoring). Aquaculture 158, 85–94. Foreman, M.G.G., 1978. Manual for Tidal Currents Analysis and Prediction. Pacific Marine Science Report 78-6. Institute of Ocean Sciences, Patricia Bay, Sidney, BC, 57 pp. Gillibrand, P.A., Gubbins, M.J., Greathead, C., Davies, I.M., 2002. Scottish Executive Locational Guidelines for Fish Farming: Predicted Levels of Nutrient Enhancement and Benthic Impact. Scottish Fisheries Research Report 63/2002. Fisheries Research Services, Aberdeen, 52 pp. Gowen, R.J., Smyth, D., Silvert, W., 1994. Modelling the spatial distribution and loading of organic fish farm waste to the
seabed. In: Hargrave, B.T. (Ed.), Modelling Benthic Impacts of Organic Enrichment from Marine Aquaculture. Canadian Technical Report on Fisheries and Aquatic Sciences, vol. 1949, pp. 19–30. Green, J.A.M., Stigebrandt, A., 2003. Statistical models and distributions of current velocities with application to the prediction of extreme events. Estuar. Coast. Shelf. 58, 601–609. Hansen, P.K., Ervik, A., Schaanning, M., Johannessen, P., Aure, J., Jahnsen, T., Stigebrandt, A., 2001. Regulating the local environmental impact of intensive marine fish farming. II. The monitoring programme of the MOM system (Modelling-Ongrowing Fish Farms-Monitoring). Aquaculture 194, 75–92. Hevia, M., Rosenthal, H., Gowen, R.J., 1996. Modelling benthic deposition under fish cages. J. Appl. Ichtyol. 12, 71–74. Krasakopoulou, E., Souvermezoglou, E., Minas, H.J., Scoullos, M., 2005. Organic matter stoichiometry based on oxygen consumption—nutrients regeneration during a stagnation period in Jabuka Pit (middle Adriatic Sea). Cont. Shelf. Res. 25, 127–142. Kloeden, P.E., Platen, E., 1992. Numerical Solution of Stochastic Differential Equations. Springer-Verlag, Berlin/Heidelberg/New York, 632 pp. Kurbanmuradov, O., Sabelfeld, K., 2000. Lagrangian stochastic models for turbulent dispersion in the atmospheric boundary layer. Bound. Lay. Meteorol. 97, 191–218. Monin, A.S., Yaglom, A.M., 1971. In: Lumley, J.L. (Ed.), Statistical Fluid Mechanics: Mechanics of Turbulence, vol. 1. The MIT Press, Cambridge, MA; London, England, p. 769. Pawlowicz, R., Beardsley, B., Lentz, S., 2002. Classical tidal harmonic analysis including error estimates in MATLAB using T TIDE. Comput. Geosci.-UK 28, 929–937. Rodean, H.C., 1996. Stochastic Lagrangian Models of Turbulent Diffusion. American Meteorological Society, Boston, 84 pp. Silvert, W., Sowles, J.W., 1996. Modelling environmental impacts of marine finfish aquaculture. J. Appl. Ichtyol. 12, 75–81. Stigebrandt, A., Aure, J., Ervik, A., Hansen, P.K., 2004. Regulating the local environmental impact of intensive marine fish farming. III. A model for estimation of the holding capacity in the Modelling-Ongrowing Fish Farm-Monitoring System. Aquaculture 234, pp. 289–261. Warner, J.C., Sherwood, C.R., Arango, H.G., Signell, R.P., 2005. Performance of four turbulence closure models implemented using a generic length scale method. Ocean Model. 8, 81–113. Wilson, J.D., 2000. Trajectory models for heavy particles in atmospheric turbulence: comparison with observations. J. Appl. Meteorol. 39, 1894–1912. Wilson, J.D., Sawford, B.L., 1996. Review of Lagrangian stochastic models for trajectories in the turbulent atmosphere. Bound. Lay. Meteorol. 78, 191–210.
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/ecolmodel
Impact of aquacultures on the marine ecosystem: Modelling benthic carbon loading over variable depth Marko Jusup ∗ , Sunˇcana Geˇcek, Tarzan Legovi´c – Boˇskovi´c Institute, Bijeniˇcka Cesta 54, 10002 Zagreb, Croatia Ruder
a r t i c l e
i n f o
a b s t r a c t
Article history:
In order to be used within Environmental Impact Assessment study, we have developed a
Received 2 August 2005
three-dimensional particle tracking model for prediction of benthic carbon loading (BCL)
Received in revised form
caused by fish farms. The model is based on stochastic differential equations for particle
15 August 2006
transport consistent with the well-known semi-empirical advection/diffusion equation. It
Accepted 24 August 2006
requires only easily obtainable input data in the form of measured current record, the source
Published on line 5 October 2006
location and a specification of local bathymetry. The model accounts for advection by longterm residual and tidal currents, turbulent diffusion, realistic bathymetry and variations in
Keywords: Fish farms
daily (monthly or yearly) emissions from fish farm. Here, we concentrate on the changes in sedimentation pattern caused by various bathy-
Benthic carbon loading
metric shapes. Examination of idealized cases reveals where and why we can expect the
Model
worst impact on benthic communities. For future reference, these results will be included
Variable depth
into guidelines for fish farming. © 2006 Elsevier B.V. All rights reserved.
1.
Introduction
Mariculture has already become very important area of food production all over the world. Potentially high profits have attracted investors in the last two decades triggering its rapid development. At the same time, many concerns have been raised about possible negative environmental impact, especially at locations where regional economy is tourism based. Main sources of environmental problems in connection with fish farming are enrichment of both water column and benthos with organic material, and physical and biological disturbances caused by cage structures (Silvert and Sowles, 1996). Therefore, in order to ensure the long-term use of farming sites Ervik et al. (1997) proposed three basic sustainability requirements and later Hansen et al. (2001) presented details of appropriate monitoring program. In these papers most concern is devoted to prevention of disappearance of benthic
∗
infauna due to excessive accumulation of particulate organic wastes. The first step toward resolution of this issue is development of a model that accurately estimates benthic nutrient loading. Hence, we present a tool, named KK3D, for prediction of dispersal of particulate organic matter from fish farms. It belongs to a broad class of Lagrangian stochastic models (LSM) which have been widely used in atmospheric sciences (see Wilson and Sawford, 1996, for a review). In oceanography the interest for them is relatively low and some of the probable reasons are discussed in Section 3. The model is based on stochastic differential equations consistent with the well-known semi-empirical advection/diffusion equation. It requires only easily obtainable input data in the form of measured current record, location of the source and the specification of local bathymetry. The model accounts for advection by long-term residual and tidal currents, turbulent diffusion, realistic bathymetry and variations in daily (monthly
Corresponding author. Tel.: +385 1 4680230. ˇ ´ E-mail addresses: [email protected] (M. Jusup), [email protected] (S. Gecek), [email protected] (T. Legovic). 0304-3800/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2006.08.007
460
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
or yearly) emissions from fish farms. Special attention is given to changes in the sedimentation pattern due to variable bathymetry. The paper is organized as follows. In Section 2 an overview of other benthic nutrient loading models is given. It is followed by a brief review of LSM theory and details of carried simulations in Section 3. Section 4 displays the results. The paper is concluded with a discussion of results in Section 5.
2.
Overview of other models
Modelling offers methods for analyzing various aspects of aquaculture like: - identification of suitable sites by GIS technique (Congleton et al., 1999), - assessment of carrying capacity for multispecies culture by coupling physical and biogeochemical models (Duarte et al., 2003), - simulation of the dynamics of an integrated experimental aquaculture plant (Buonomo et al., 2005), long before the production is initiated. Environmental impact estimation of marine fish farms, on both local and regional level, is no exception. So far, much attention was given to particle tracking simulations that predict distribution of particulate organic wastes on the seabed (i.e. the sedimentation pattern or the deposition footprint). One of the simpler approaches (Gowen et al., 1994) is based on the observation that the dispersion distance D is related to current speed U, water depth H and settling velocity of pellets w by D=
UH . w
(1)
It represents the distance that particles travel from the cage before settling on the bottom. The presence of water depth in Eq. (1) makes the “Gowen” method suitable to include variable bathymetry. However, it was Hevia et al. (1996) who first extended the model to account for bottom topography and vertical variations in current speed. The importance of turbulent diffusion in horizontal dispersion was also realized, resulting in an additional term in the relation for dispersion distance (Gillibrand et al., 2002) UH D= + w
2KH . w
where K is the turbulent diffusion coefficient. Stigebrandt et al. (2004) adopted somewhat different approach. They estimate the current speed variance 2 from the measurements and use it to determine the dispersion length T = H/w of a given site. For each dispersion length T there is a pre-calculated, dimensionless, normalized loading (sedimentation) function = (r) which relates emission from the net pens F1 with benthic carbon loading F2 through the relationship F2 (r) = (r)F1 . The distance from the cage center is denoted by r. Note how this procedure should be computationally inexpensive once
the loading functions have been determined. However, that would be an impossible task without the fact that current velocities appear to be approximately normally distributed (Green and Stigebrandt, 2003). DEPOMOD, developed by Cromey et al. (2002a), is a tool based on direct particle tracking. It was validated in sediment trap experiments and showed good agreement with field data with accuracy of ±20 and ±13% for a dispersive and depositional site, respectively. The model successfully links solids accumulation to benthic indices and total individual abundance. Two processes determine the path of a particle: advection by measured currents and turbulent diffusion modelled with random walk. This approach is analogous to the methods used in KK3D. DEPOMOD is accompanied by the validated resuspension model (Cromey et al., 2002b) which is a rare occurrence within the field. All the mentioned particle tracking models are characterized by physical simplicity and reasonable predictive capability. However, none of them is equivalent to the advection/diffusion equation which governs the tracer dispersion from a given source. The use of LSM enabled us to achieve that goal. Consequently, KK3D is a tool based on proven physical principles and underlying stochastic differential equations are quite simple to implement as a numerical simulation on a PC. Furthermore, recent inclusion of the adsorption term naturally accounted for nutrient loss caused by leaching. Information on the chemical composition of organic matter in faeces reaching the bottom allows the calculation of oxygen demand (g O2 /g C). For that purpose stoichiometric analysis similar to Krasakopoulou et al. (2005) was adopted. Thus, KK3D addresses the chemical and biological processes which determine the final fate of nutrients emitted from a fish farm. However, full description of these additions is beyond the scope of this paper. We cannot claim that KK3D resolves all questions regarding the sedimentation pattern. The model user should be aware of an important limitation because there is no resuspension of sediment. It is a physical process which, beside the mentioned chemical and biological processes, affects the final distribution of organic waste on the seabed. Resuspension is usually split into four compartments: erosion, transport, deposition and consolidation. Erosion occurs when the shear stress on the bottom = U∗2 becomes greater than the critical erosion shear stress ce . Here, is the density of seawater and U* is the bed shear velocity which is estimated from the measured current speed near the bottom Uz assuming a logarithmic profile U∗ =
Uz . ln(z/z0 )
Thus, the hydraulic roughness length z0 , as an important characteristic of the sediment type, plays a role in the formation of sedimentation pattern. Transport is modelled by the same equations used for determining the initial distribution of waste material. Deposition occurs when the shear stress on the bottom becomes lower than the critical deposition shear stress cd . It is assumed that cd < ce so that erosion and deposition are mutually exclusive. Consolidation is a function of
461
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
time spent on the seabed. The longer the particle stays on the bottom, probability of being eroded approaches zero. To conclude, only a full resuspension model can account for properties of particular sediment type and this has not yet been included into the KK3D tool. Finally, we fell that the oxygen demand should be related to the flux of oxygen into the sediment and oxygen consumption by the microbial degradation. A promising attempt to resolve that problem was made by Stigebrandt et al. (2004).
3.
Methods
It is a common practice to breakdown the sea current velocity into two parts: (1) slowly varying component due to tidal, thermohaline and wind forcing and (2) random component related to turbulence that quickly changes in time. Former part is responsible for displacement of waste material from the source (advection) while the latter part causes the dispersion (turbulent diffusion). Therefore, it is tempting to adopt the semi-empirical approach (Monin and Yaglom, 1971) and model the transport of any passive tracer emitted from an aquaculture by using the advection/diffusion equation: ∂C ∂C ∂ = + vi (x, t) ∂t ∂xi ∂xi
Kij (x, t)
∂C ∂xj
,
vi (x, t) + vz ıi3 +
∂Kij (x, t) ∂xj
dt +
vi C = Kij
∂C . ∂vj
Particle tracking model described above requires two types of input data:
(2)
where C is a tracer concentration, vi , i = 1, 2, 3 the current velocities and Kij , i, j = 1, 2, 3 are the eddy diffusivity coefficients. This is an Eulerian approach and we can replace it by a Lagrangian technique called Random Displacement Model (RDM) which is consistent with Eq. (2) (Rodean, 1996). However, in order to predict the distribution of particulate organic material on the seabed we must calculate the paths of heavy particles and not the passive tracer. According to Wilson (2000), RDM easily resolves this problem by superimposing to it the particle settling velocity vz . Thus, the position of a particle in RDM is given by stochastic differential equations:
dxi =
2005). Furthermore, there is a problem of so-called trajectory crossing. Initially the heavy particle position will coincide with the position of some fluid element. As the time passes, particle drifts away from the fluid element because of its settling rate, leaving the eddy in which correlation with initial velocity still exists. Therefore, it is suggested to reduce the heavy particle correlation timescale with respect to the fluid element timescale (e.g. Wilson, 2000), but this is a heuristic approach. There is also a possibility that due to finite inertia, heavy particle is unable to follow rapid velocity fluctuations, resulting in a different velocity variance compared to the fluid element. Our objective was to construct a model that is usable in many different locations prior to installation of an aquaculture. In such situations, at best, measured current profiles are available as input data. So to avoid unnecessary complexities and account for dispersion, one must impose reasonable assumptions about joint moments of concentration and velocity fluctuations. The simplest and widely used example is the eddy diffusion hypothesis
2Kij (x, t)dWj ,
(3)
where dWi , i = 1, 2, 3 are the standard Wiener processes. Sedimentation pattern is expressed in form of benthic carbon loading (BCL) in gC/m2 day. Wilson (2000) classifies RDM as a “zeroth-order” LSM noting its main weakness; if xi (t) is treated as Markovian random variable, the model cannot be valid close to the source, that is for travel times shorter than the velocity correlation timescale. The more advanced, “first-order” LSM, called Langevin Equation Models (LEM), have been used extensively in atmospheric research owing much of their popularity to analytical description of turbulence properties in the surface boundary layer. For example, Kurbanmuradov and Sabelfeld (2000) have shown how a number of simplifying assumptions lead to unique selection of three-dimensional LEM for turbulent diffusion in the atmosphere. This turns out to be an impossible task in oceanographic application (Brickman and Smith, 2002) where turbulent quantities are provided by circulation models with incorporated appropriate closure scheme (e.g. Warner et al.,
1. Measured current record from a representative location must be obtained. Slowly varying tidal components as well as long-term residual mean should be extracted from such a record while short-period fluctuations related to turbulence are of no interest here since we adopted the eddy diffusion assumption. Processing of current meter data by classical harmonic analysis (Pawlowicz et al., 2002) provides satisfactory results returning the following information about each tidal component: - is the frequency of tidal component in cycles per hour (cph), - a is the major semi-axis of tidal current ellipse (cm/s), - b is the minor semi-axis of tidal current ellipse (cm/s), - is the inclination of tidal current ellipse, - g is the Greenwich phase, - V(t0 ) is the astronomical argument at the moment t0 , - u is the nodal phase modulation, - f is the nodal amplitude correction. From these parameters, the model generates tidal current according to following relations (Foreman, 1978): a+ =
a+b f, 2
a− =
a−b f 2
ε+ = − g + (u + V(t0 )),
ε− = + g − (u + V(t0 ))
(4a) (4b)
v1 (t) = a+ cos(ε+ + 2 t) + a− cos(ε− − 2 t)
(4c)
v2 (t) = a+ sin(ε+ + 2 t) + a− sin(ε− − 2 t).
(4d)
Thanks to Eqs. (4a)–(4d) KK3D is capable of truly forecasting the sedimentation pattern and not only hind-casting as it is the case when measured current record is directly loaded into
462
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
the model. Extreme events, like storms, can be modelled by modifying the residual current for a short period of time. Modification should be made according to measurements or results from an appropriate hydrodynamic model. Also, activities like net cleaning are easily included into the prediction. The user should simply increase farm output during the desired time period. 2. Rectangular grid containing bottom topography is also needed. Grid cell resolution is arbitrary and depends only on available data from a location of interest. Cromey et al. (2002a) used the resolution of 10 or 25 m depending on the size of the sedimentation pattern and in realistic application one can hardly expect higher precision due to instrument limitations. Here, we have concentrated on emissions from a single net pen, which cover relatively small area on the seabed. Therefore, slightly denser grid is used for obtaining the results discussed in this paper. To start particle tracking we also need the magnitude of settling velocity vz . The velocity is chosen randomly for each particle from the normal distribution with the mean = 3.2 cm/s and the standard deviation = 1.1 cm/s (Cromey et al., 2002a). Approximate solutions to Eq. (3) are obtained by means of the Euler–Maruyama method or, as Kloeden and Platen (1992) put it, an Ito–Taylor scheme of strong order 0.5:
xi (t + t)xi (t) +
vi (x, t) + vz ıi3 +
∂Kij (x, t) ∂xj
t +
2Kij (x, t) t j .
Here, t is the time-step and j ∈ N(0, 1). In the atmospheric surface boundary layer, it is convenient to choose coordinate axes so that the x1 -axis coincides with the direction of the mean wind velocity. Then, the crosswind velocity fluctuations are symmetrical with respect to the planes x1 = 0 and x3 = 0, and we have v1 v2 = 0 and v2 v3 = 0. In addition, since the Ox1 , Ox2 and Ox3 are clearly the preferred axes of the flow they are often considered as principal axes of the diffusivity tensor. Within the sea, the direction of mean flow can change with depth making the above reasoning inapplicable. However, Monin and Yaglom (1971) note that “all qualitative predictions of semi-empirical diffusion theory can hardly depend on the presence or absence of terms with coefficients K13 and K31 ” (authors considered the inclusion of these terms into Eq. (2)). Accordingly, our model includes only diagonal elements of the diffusivity tensor. Besides, these coefficients are depth dependent to account for the existence of surface and bottom boundary layers. Their mean values correspond to the values observed in Scottish inshore waters (K11 = K22 = 0.1 m2 /s and K33 = 0.001 m2 /s) and are used for modeling by Gillibrand et al. (2002) and Cromey et al. (2002a). Thus, based on published results, we expect that the key properties of sedimentation pattern such as its shape and size, the location of maximum BCL, etc., can be estimated with acceptable accuracy. Initial particle position in x1 Ox2 plane is chosen randomly from a uniform distribution assuming that the fish can be located anywhere within the cage limits immediately prior to defecation. In all simulations the cage radius was set to 20 m and the center was placed at x1 = x2 = 0. Let us note that fix-
Table 1 – Tidal ellipse parameters for constituents used in all simulations Tidal constituent
O1 K1 M2 S2
Frequency (cph) 0.0387307 0.0417807 0.0805114 0.0833333
Major semi-axis (cm/s) 2.00 5.32 8.50 5.20
Minor semi-axis (cm/s) −0.23 −0.36 −0.32 −0.27
Negative minor semi-axis indicates clockwise rotation of tidal current phasor. Inclination was set to zero. Note the highly polarized tidal currents typical for the Adriatic Sea.
ing the starting position (for example, to the cage center) does not affect the sedimentation pattern at all in sufficiently deep water. In such conditions, the deposition area is much larger than the cage area and the exact initial position of the particle is hardly relevant. However, in shallow water fixing the starting position can result with under-prediction of sedimentation pattern size. The model is far more sensitive to the choice of initial vertical coordinate x3 because dispersion is depth dependent, i.e. deposition area increases rather quickly with depth. Since there is observational evidence that defecation occurs near the surface immediately after feeding, we have chosen x3 = −5 m for all particles. The number of tracked particles in each simulation was set to 2 × 105 since no changes in results occurred with higher values. Output data are passed to visualization module, which produces contour plots revealing details of predicted sedimentation pattern. Unless we simulate exceedingly large number of sample paths, the stochastic nature of employed algorithm will cause appearance of small patches in final images. Clearly, these patches are insignificant as they disappear with large number of tracked particles. To avoid spending too much computational time we recommend the removal of model patchiness by a two-dimensional wavelet analysis. In this study, dispersion is simulated over 1-day period, but that value can easily be changed in order to meet specific requirements. For example, sometimes it is useful to predict dispersion when strong currents towards the shore dominate the aquaculture area. The model also accounts for uneven emission rate of waste material during the simulated period. This is necessary as we expect heavier impact in direction of prevailing current when emission rate reaches maximum. Here, we have concentrated on changes in sedimentation pattern induced by various bathymetric features. Consequently, for all examined cases the residual and tidal currents are the same. In Table 1 information on used tidal constituents is given. Values are typical for the Adriatic Sea. One should keep in mind that all studied cases are idealized.
4.
Results
In order to investigate the validity of described model we have run a series of simulations that were supposed to satisfy the following criteria: (i) in case of disabled advection, the model should reproduce the results of Stigebrandt et al. (2004) illustrated in dispersion submodel section of their paper and (ii) in
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
463
Fig. 1 – Schematic view of bathymetry over a transect y = 0 for all simulations: (a) flat seabed, (b) 2.5◦ slope, (c) 5◦ slope, (d) 5 m deep pit centered at x = −100 m and (e) 5 m high hill centered at x = −60 m. Arrows in (a) indicate direction and magnitude of residual current.
case of bottom slope, the model should reproduce the effects reported by Hevia et al. (1996). Both criteria were successfully fulfilled; instance (i) is not further discussed here because it is not directly related to the main topic, but (ii) is explained in detail below. Computer simulation results are illustrated in Figs. 1 and 2. Fig. 1 shows the crosssectional view of bathymetry along the line y = 0 as well as the position of net pen for all studied cases. Contour plots in Fig. 2 display the corresponding sedimentation patterns. The first simulation was made for the flat seabed and it is considered a reference point to which all other cases are compared. Before examining the changes in sedimentation pattern caused by several different bottom topographies, we must introduce few important concepts. Let us denote the maximum BCL by m so that the values of BCL fall within the segment S = [0, m]. Subdivision of S into six subintervals S1 = [0, x1 ], S2 = [x1 , x2 ], . . ., S6 = [x5 , m] (where S1 , . . ., S5 are equally long) provides the basis for contour plots shown in Fig. 2, that is, isolines (lines of constant BCL) are drawn at levels xi (i = 1, . . ., 5). Note how such procedure exposes
the details of each sedimentation pattern and allows one to define the maximum deposition region as a set of coordinates M = {(x, y) : BCL(x, y) ∈ S6 }. If surface z = z(x, y) represents the seabed (i.e. z is the depth at (x, y)), it is possible to approximate an area A corresponding to the region M. Parameter A, named the maximum deposition area, provides valuable information for comparing various test cases. Other features revealed by contour plots are the shape of M and the position where maximum BCL occurs, i.e. the coordinates (x0 , y0 ) such that m = BCL(x0 , y0 ). The impact on the bottom is considered negligible if the BCL level is below 0.1 gC/m2 day. Table 2 is a convenient summary for the presentation of results that follow below. For the reference simulation, constant depth level was set to 35 m (Fig. 1a). Sedimentation pattern is elongated in agreement with the polarized tidal current and displaced in response to residual current (Fig. 2a). Maximum deposition region is split into two parts what is, again, a consequence of polarized tidal current which reverses its direction throughout a day. With assumed emission of 10 kgC/day in the form of
464
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
Fig. 2 – BCL (in gC/m2 day) for all simulated test cases: (a) flat seabed, (b) 2.5◦ slope, (c) 5◦ slope, (d) 5 m deep pit centered at x = −100 m and (e) 5 m high hill centered at x = −60 m. Maximum BCL location is marked with ‘x’. Dashed circle is a projection of cage on the bottom.
feces, the maximum BCL reaches the value of 1.17 gC/m2 day at the distance of 84 m from the cage center. The total area affected by organic waste is 21,840 m2 and the area of maximum deposition equals 1704 m2 .
Even a small variation in bathymetry can cause major changes in the sedimentation pattern. This is shown in the simulation for 2.5◦ bottom slope (Fig. 1b). Note that this is the test case from the criterion (ii) for the validity of our model
Table 2 – A summary of the simulation results Maximum BCL (gC/m2 day) 1.17 1.55 1.92 1.20 1.56
Maximum BCL position (m) −84, 4 −72, 0 −66, −4 −114, 0 −54, 0
Total affected area (m2 )
Maximum deposition area (m2 )
21840 20828 20082 22165 21925
1704 2474 2722 1325 1810
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
as mentioned above. The depth range caused by a slope is 22–48 m with the mean of 35 m. Sedimentation pattern in Fig. 2b discloses the restricted deposition towards the shallow water because the tilt blocks the particles from reaching the coastline. The value of maximum BCL increases to 1.55 gC/m2 day and is attained at 72 m from the cage center. Total area influenced by organic material decreases to 20,828 m2 , but the maximum deposition area grows up to 2474 m2 . Doubling the bottom slope to 5◦ makes the changes in sedimentation pattern more prominent (Fig. 1c). Again, the mean depth is 35 m, but steeper slope causes wider depth range over the simulation area (9–61 m). At 150 m from the cage center towards the shallow area there is no noticeable deposition (Fig. 2c). For the reference simulation and in case of gentler slope the same happens at 210 and 174 m from the cage center, respectively. Hence, the deposition in shallow water is more restricted by higher slope steepness. The maximum BCL reaches the value 1.92 gC/m2 day and is closer to the net pen than previously (x0 = −66 m). The total area disturbed by organic waste reduces to 20,082 m2 and the maximum deposition area is 2722 m2 . Following bottom topographies were designed to illustrate how by neglecting relatively small bathymetric features one is lead to miss-prediction of all important parameters. More precisely, the model can produce wrong shape of the maximum deposition region (M) and give incorrect values for its size (A), the maximum BCL (m) and the corresponding coordinates (x0 , y0 ). Consider the 5 m deep pit with inner wall steepness of less than 10◦ (Fig. 1d). It is located on an otherwise flat bottom that is 35 m deep (as in the reference simulation). The pit is centered at 100 m from the cage center, meaning that it does not coincide with original maximum deposition region. However, the presence of such pit has important consequences on sedimentation pattern and most notable is the displacement (of one part) of M further away from the net pen (Fig. 2d). The maximum BCL of almost 1.20 gC/m2 day is now located at x0 = −114 m. At the location where maximum BCL occurs originally (x0 = −84 m), the predicted value has changed to 0.68 gC/m2 day; not even close to the highest loading. Seabed area that is worst affected is somewhat reduced (A = 1325 m2 ) as a consequence of better dispersal over the pit. Exactly the opposite changes are induced by 5 m high hill centered at the distance of 60 m from the net pen (Fig. 1e). The hill, as a projection from the bottom, prevents particles from reaching their original settling location. For this reason (a part of) the maximum deposition region is now displaced closer to the cage (Fig. 2e). Also, being the highest point on the bottom, the hill makes the dispersion less effective. As a result, the maximum BCL increases to 1.56 gC/m2 day and that value is recorded at x0 = −54 m; near the hill top. At the location where the maximum BCL occurs originally (x0 = −84 m), the predicted value has changed to 0.74 gC/m2 day. Thus, ignoring the hill we would completely miss both parameters: the highest BCL m and corresponding coordinates (x0 , y0 ). Finally, there are no major changes in the total affected area (21,925 m2 ) and maximum deposition area (1810 m2 ).
5.
465
Discussion
We have presented the capabilities of a model for prediction of BCL caused by fish farms. The model is based on strong, yet relatively simple, physical grounds and can account for number of features important in shaping the organic waste sedimentation pattern. Among these features, special attention is given to modelling the dispersion over variable bathymetry. For this purpose we have compared the results from several simulations with the sedimentation pattern obtained for the flat seabed, holding other environmental variables constant. The inevitable conclusion is that approximating the realistic bathymetry with the uniform bottom can result in serious error. It is possible to predict wrong shape of the maximum deposition region, miscalculate its area, give incorrect value and location of the maximum BCL and even fail to estimate the total area exposed to particulate waste. Thus, as others have stated before, necessity for more robust models exists, but the balance between complexity and available data from field must be retained. Accordingly, the presented model is based on data that can easily be recorded in practice and (even more importantly) with limited expenses. Modular programming concept, popular because it easily accommodates further extensions of the model (Cromey et al., 2002a; Stigebrandt et al., 2004), has been applied here in order to finally produce user-friendly application for fish farm installation and management. So far, the described BCL module has been used in several Environmental Impact Assessment studies (required by the Croatian Ministry of Environmental Protection, Physical Planning and Construction) for farms with production that exceed 50 t/year. In future we expect the model to become the standard, not only in the selection of suitable locations, but also for the improvement of farming practices. This implies achieving the maximum production that ensures long-term use of a site. For example, one way of approaching that goal is linking BCL to oxygen consumption in order to predict the change from aerobic to anaerobic conditions in the sediment. As this change may affect the farm performance, it would be useful to predetermine the optimum tonnage for the net pens (i.e. carrying capacity) at which production levels do not decrease. Future development may proceed in two directions: (1) refinement of physical aspects of the model and (2) conversion of predicted quantities to biochemical parameters. Regarding point (1) the most obvious step is constructing the module for post-depositional changes like the slumping of material down the steep slopes or resuspension of sediment into the water column. Another possibility is coupling of dispersion predictions with full hydrodynamic model that can resolve currents horizontally and calculate eddy diffusivities. However, practical use of such model is questionable if we consider the input data necessary for running it. Perhaps the most important issue of point (2) is connection of BCL with oxygen consumption in sediments as explained above. Local and regional water quality is also an important parameter. On the local level the water should have sufficient oxygen concentration and low concentrations of potentially harmful substances. Regional water quality implies insignificant contribution of mariculture to algal blooms which occur as a consequence of nutrient
466
e c o l o g i c a l m o d e l l i n g 2 0 0 ( 2 0 0 7 ) 459–466
enrichment. Some of these topics have already been addressed in the literature, but definitive solutions are yet to be reached.
Acknowledgments This research has been supported by the Croatian Ministry for Science, Education and Sport, the FP6 Project on ecosystem approach to sustainable aquaculture and the NorwegianCroatian Project on the particularly sensitive sea area.
references
Brickman, D., Smith, P.C., 2002. Lagrangian stochastic modelling in coastal oceanography. J. Atmos. Ocean. Tech. 19, 83–99. Buonomo, B., Falcucci, M., Hull, V., Rionero, S., 2005. A mathematical model for an integrated experimental aquaculture plant. Ecol. Model. 183, 11–28. Congleton Jr., W.R., Pearce, B.R., Parker, M.R., Beal, B.F., 1999. Mariculture siting: a GIS description of intertidal areas. Ecol. Model. 116, 63–75. Cromey, C.J., Nickell, T.D., Black, K.D., 2002a. DEPOMOD-modelling the deposition and biological effects of waste solids from marine cage farms. Aquaculture 214, 211–239. Cromey, C.J., Nickell, T.D., Black, K.D., Provost, P.G., Griffiths, C.R., 2002b. Validation of a fish farm waste resuspension model by use of a particulate tracer discharged from a point source in a coastal environment. Estuaries 25, 916–929. Duarte, P., Meneses, R., Hawkins, A.J.S., Zhu, M., Fang, J., Grant, J., 2003. Mathematical modelling to assess the carrying capacity for multi-species culture within coastal waters. Ecol. Model. 168, 109–143. Ervik, A., Hansen, P.K., Aure, J., Stigebrandt, A., Johannessen, P., Jahnsen, T., 1997. Regulating the local environmental impact of intensive marine fish farming. I. The concept of the MOM system (Modelling-Ongrowing Fish Farms-Monitoring). Aquaculture 158, 85–94. Foreman, M.G.G., 1978. Manual for Tidal Currents Analysis and Prediction. Pacific Marine Science Report 78-6. Institute of Ocean Sciences, Patricia Bay, Sidney, BC, 57 pp. Gillibrand, P.A., Gubbins, M.J., Greathead, C., Davies, I.M., 2002. Scottish Executive Locational Guidelines for Fish Farming: Predicted Levels of Nutrient Enhancement and Benthic Impact. Scottish Fisheries Research Report 63/2002. Fisheries Research Services, Aberdeen, 52 pp. Gowen, R.J., Smyth, D., Silvert, W., 1994. Modelling the spatial distribution and loading of organic fish farm waste to the
seabed. In: Hargrave, B.T. (Ed.), Modelling Benthic Impacts of Organic Enrichment from Marine Aquaculture. Canadian Technical Report on Fisheries and Aquatic Sciences, vol. 1949, pp. 19–30. Green, J.A.M., Stigebrandt, A., 2003. Statistical models and distributions of current velocities with application to the prediction of extreme events. Estuar. Coast. Shelf. 58, 601–609. Hansen, P.K., Ervik, A., Schaanning, M., Johannessen, P., Aure, J., Jahnsen, T., Stigebrandt, A., 2001. Regulating the local environmental impact of intensive marine fish farming. II. The monitoring programme of the MOM system (Modelling-Ongrowing Fish Farms-Monitoring). Aquaculture 194, 75–92. Hevia, M., Rosenthal, H., Gowen, R.J., 1996. Modelling benthic deposition under fish cages. J. Appl. Ichtyol. 12, 71–74. Krasakopoulou, E., Souvermezoglou, E., Minas, H.J., Scoullos, M., 2005. Organic matter stoichiometry based on oxygen consumption—nutrients regeneration during a stagnation period in Jabuka Pit (middle Adriatic Sea). Cont. Shelf. Res. 25, 127–142. Kloeden, P.E., Platen, E., 1992. Numerical Solution of Stochastic Differential Equations. Springer-Verlag, Berlin/Heidelberg/New York, 632 pp. Kurbanmuradov, O., Sabelfeld, K., 2000. Lagrangian stochastic models for turbulent dispersion in the atmospheric boundary layer. Bound. Lay. Meteorol. 97, 191–218. Monin, A.S., Yaglom, A.M., 1971. In: Lumley, J.L. (Ed.), Statistical Fluid Mechanics: Mechanics of Turbulence, vol. 1. The MIT Press, Cambridge, MA; London, England, p. 769. Pawlowicz, R., Beardsley, B., Lentz, S., 2002. Classical tidal harmonic analysis including error estimates in MATLAB using T TIDE. Comput. Geosci.-UK 28, 929–937. Rodean, H.C., 1996. Stochastic Lagrangian Models of Turbulent Diffusion. American Meteorological Society, Boston, 84 pp. Silvert, W., Sowles, J.W., 1996. Modelling environmental impacts of marine finfish aquaculture. J. Appl. Ichtyol. 12, 75–81. Stigebrandt, A., Aure, J., Ervik, A., Hansen, P.K., 2004. Regulating the local environmental impact of intensive marine fish farming. III. A model for estimation of the holding capacity in the Modelling-Ongrowing Fish Farm-Monitoring System. Aquaculture 234, pp. 289–261. Warner, J.C., Sherwood, C.R., Arango, H.G., Signell, R.P., 2005. Performance of four turbulence closure models implemented using a generic length scale method. Ocean Model. 8, 81–113. Wilson, J.D., 2000. Trajectory models for heavy particles in atmospheric turbulence: comparison with observations. J. Appl. Meteorol. 39, 1894–1912. Wilson, J.D., Sawford, B.L., 1996. Review of Lagrangian stochastic models for trajectories in the turbulent atmosphere. Bound. Lay. Meteorol. 78, 191–210.