Ecosystem modelling for ecosystem-based management of bivalve aquaculture sites in data-poor environments

Although models of carrying capacity have been around for some time, their use in aquaculture management has been limited. This is partially due to the cost involved in generating and testing the models. However, the use of more generic and flexible models could facilitate the implementation of modelling in management. We have built a generic core for coupling biogeochemical and hydrodynamic models using Simile (www.simulistics.com), a visual simulation en vironment software that is well-suited to accommodate fully spatial models. Specifically, Simile integrates PEST (model-independent parameter estimation, Watermark Numerical Computing, www. pesthomepage.org), an optimization tool that uses the Gauss-Marquardt-Levenberg algorithm and can be used to estimate the value of a parameter, or set of parameters, in order to minimize the discrepancies between the model results and a dataset chosen by the user. The other critical aspect of modelling exercises is the large amount of data necessary to set up, tune and ground truth the ecosystem model. However, ecoinformatics and improvements in remote sensing procedures have facilitated acquisition of these datasets, even in data-poor environments. In this paper we describe the required datasets and stages of model development necessary to build a biogeochemical model that can be used as a decision-making tool for bivalve aquaculture management in data-poor environments.


INTRODUCTION
Coastal zones provide two-thirds of global ecosystem services (Costanza et al. 1997), but anthropogenic stresses such as over-fishing, pollution and the direct and indirect impacts of climate change compromise their sustainability (Hughes et al. 2005). Maintenance of ecosystem functioning to provide the services humans require is the goal of ecosystembased management (EBM; McLeod et al. 2005). EBM necessitates the incorporation of ecosystem resili-ence, i.e. the capacity of a system to absorb disturbance and reorganize while undergoing change, so as to retain essentially the same functions, structure, identity and feedbacks (Walker et al. 2004). These concepts have been applied to coastal aquaculture, where there are concerns that feeding and wastecultured animals such as bivalves or fish alter particle and nutrient fluxes to the detriment of benthic health. In the context of coastal bivalve aquaculture, the filtration activity removes phytoplankton and has the potential to impact top-down control of primary pro-duction, which in certain locations could exert a positive effect by mitigating eutrophication (Coen et al. 2007).
These effects have been studied by means of carrying capacity, which is defined by 4 components: physical, ecological, production and social (for production carrying capacity: McKindsey et al. 2006, Ferreira et al. 2012. Grant et al. (2007) and  have combined the concepts of ecosystem functioning and resilience in an alternative definition of carrying capacity, i.e. the bivalve stocking density at which growth is not food limited and/or some measure of ecosystem health is not compromised, which is applied in this study using ecosystem modelling. Ecosystem models are powerful tools for exploring carrying capacity, allowing the study of stocks, energy fluxes and potential interactions in complex ecosystems (Dowd 2005). Models integrate time and space, which is critical for understanding ecological dynamics and therefore how natural systems provide ecosystem services (Palmer et al. 2004). In addition, scenario building allows the exploration of future situations where unanticipated stressors generate new risks or opportunities and is thus an important tool for managing those changes (Nobre et al. 2010).
Consideration of resilience may be regarded as a perspective for the analysis of social-ecological systems wherein optimization is an outcome-oriented tool used to make rational and transparent decisions about a well-defined problem (Fischer et al. 2009). In ecosystem modelling, optimization works to make a system as efficient as possible. However, common to all applications of optimization in conservation ecology is the quantitative identification of a conservation problem, that is, the precise definition of the limits at which ecosystem health is not compromised , Fischer et al. 2009). These limits are known as tipping points (Fig. 1), or the critical thresholds at which a small perturbation can qualita-tively alter the state or development of a system (Lenton et al. 2008). When a perturbation goes be yond a tipping point the resilience is exceeded and the system reorganizes (Crowder & Norse 2008), altering ecosystem functioning and consequently ecosystem services. Therefore, precautionary definitions of these tipping points are crucial in order to optimize ecosystem functioning. In this way, EBM has been defined as a decision-making process based on alteration of ecosystem functioning within the bounds of natural variation . This criterion has been applied to aquaculture sites in order to retrospectively compare density-dependent scenarios and carrying-capacity variation in different years (Filgueira & Grant 2009) and to prospectively predict the optimal mussel standing stock biomass that would satisfy a carrying-capacity criterion in a hypothetical aquaculture scenario (R. Filgueira unpubl. data).
Together, ecosystem modelling and optimization are the ideal tools for marine EBM, because they allow exploration of resilience and tipping points as well as efficiency in providing ecosystem services without compromising the sustainability of the ecosystem. In the context of mussel aquaculture, carrying-capacity models have been around for many years but their use in management has been limited. This is partially due to the cost involved in creating and testing the model, which includes development of the ecosystem model itself and collecting input and validation data. Increasing the model complexity attempts to bolster ecological realism, but imperfect knowledge of relationships and parameters may also lead to greater scientific uncertainty (FAO 2008). Therefore, the assumption that extra details are beneficial appears to be flawed when applied at the scale and number of dimensions of end-to-end models (Fulton 2010). This implies that modelling should re strict its focus to relevant components and critical dynamics, which must be defined based on the management question to be addressed, available data, the important system features (including forcing conditions) and the appropriate scales (FAO 2008, Fulton 2010. Bi valve carrying-capacity models should in clude both ecosystem-and local-scale information in order to consider transport processes and submodels for organism and population levels (Smaal et al. 1997  term processes in shallow coastal ecosystems. The bio geochemical model is constructed in Simile (www. simulistics.com), a visual modelling en vironment (VME) that includes an optimization tool (PEST, model-independent parameter estimation, Watermark Numerical Computing, www. pesthomepage. org). In Simile, the user can define the components and processes by means of a graphical interface and use PEST to optimize the parameters of the model. In addition to the flexibility in defining and constructing the biogeochemical model, the coupled scheme has the flexibility to accept different scales in time and space. The temporal scale can be easily modified during the definition and construction of the model. Regarding spatial scales, the coupled model can be easily designed as a multiple box or a fully spatial model, depending on the purpose of the modelling exercise. Given the large influence of hydrodynamics in marine systems, detailed spatial resolution is often desired, for example, when a geographic feature exerts a large influence on circulation.
Although the generic trophic structure of the model can be re-used at different sites, circulation models are site-specific and must be tailored for each new application, as are the appropriate datasets necessary to set up, tune and groundtruth the ecosystem components. Information used to characterize and initialize the domain, far-field time series used to force the model and datasets used to groundtruth the results are critical for its execution and validation. Despite the value of model predictions and scenario building, it is impractical to design a new research program for every application of carrying-capacity models, especially given the growth of shellfish farming. These tasks can however be facilitated by novel tools and techniques. For example, ecoinformatics aims to facilitate environmental research and management by maintaining access to databases which are often established by government agencies. Similarly, satellite remote sensing of ocean color is a powerful tool to supplement field data in providing time series for forcing and groundtruthing models. The use of satellite images for the study of the ocean is continuously providing better spatial resolution and new products that allow the creation of time series of, e.g., temperature, salinity, chlorophyll, primary production, turbidity and dissolved organic matter. However, satellite remote sensing near the coast is subject to limitations due to the optical complexity of coastal waters (Moses et al. 2009). Recently, satellite remote sensing has been used to feed individualbased models in order to predict bivalve growth (Thomas et al. 2011, Filgueira et al. 2013).
Based on this scenario, we used the available physical− biogeochemical coupling scheme published by Filgueira et al. (2012) to construct a fully spatial ecosystem model based on PNZ-type dynamics (Phytoplankton− Nutrients−Zooplankton), with the addition of mussel and seston submodels. The model has been simplified to consider the most important components and processes; this significantly reduces the cost of the study in terms of model complexity and required datasets. We have sought to streamline the process of model implementation for data-poor environments, which we define as sites with limited environmental data but a reliable description of the aquaculture activity (i.e. culture practice and bivalve growth), which is required for model set up and validation. Data limitations on the ambient environment are overcome through the use of remote sensing, hydrologic modelling of the watershed and ecoinformatics from government databases. The technical aspects of data collection are extensively described in the 'Materials and methods' section for a general scenario, so that they can be used as guidelines for further studies in aquaculture research. We demonstrated in a case study of St. Ann's Harbour (eastern Canada) that these tools are sufficiently developed to make up for a lack of field data and create a model of carrying capacity for bivalve aquaculture. Moreover, we showed that ecosystem resilience can be ex plored in such a model by means of different 'what-if' simulations using variable mussel stocking densities. Finally, we used the model to optimize sustainable mussel production as defined by conservation of ecosystem function, em ploying natural variation in chlorophyll as a benchmark. These prerequisites were used to: (1) Review and describe novel techniques such as satellite remote sensing, geographical information system (GIS) tools and optimization to collect and use data for developing ecosystem models in aquaculture data-poor sites.
(2) Apply these techniques to a case study of St. Ann's Harbour (eastern Canada) as proof of concept in order to evaluate current aquaculture activity from a sustainable standpoint.

MATERIALS AND METHODS
This section is organized in 2 sections: one section describes a general scenario in which aspects related to data collection, physical−biological coupling and optimization tools are described and a second section follows the same structure describing the case study of St. Ann's Harbour.

Data collection for modelling bivalve aquaculture sites
The information required to initialize and run a model can be grouped into 3 categories: far-field, model-domain and groundtruthing datasets. The following datasets (see Table 1 for a summary) are among the minimum required to determine the carrying capacity of a bay with mussel aquaculture activity.

Far-field
The datasets in the far field yield the characteristics of the environment beyond the model domain, and therefore are not affected by the domain itself. Conceptually, the effects of processes that occur in the model domain are rapidly diluted and would not be expected to influence the far field. These boundary conditions are used to force the model via time series covering the simulation period as follows: • Dissolved nutrients are the inorganic compounds that support phytoplankton populations. In general, pelagic estuarine productivity tends to be limited by nitrogen (Boynton & Kemp 2008); therefore, nitrogen time series in open waters are required for the model. In some ecosystems, phosphorous and silicate could also be limiting; therefore, in these specific cases, characterization of these nutrients is also necessary. The model domain in coastal environments may also be connected to rivers. This is crucial for nutrient input given that rivers are potentially rich sources of nutrients, especially in agricultural areas. Time series of nutrients can be generated either by regional climatology or land-used modelling.
• Temperature exerts a significant effect on several biological processes; for example, time series are required to adjust physiological rates to environmental conditions. These temperature records can be obtained by satellite remote sensing (Minnett et al. 2004).
• Chlorophyll concentration is used as a proxy for phytoplankton abundance. Phytoplankton are the base of the trophic web and the primary source of food for bivalves. Time series of chlorophyll content can be generated by satellite remote sensing (O'Reilly et al. 1998). As a component of model simplicity, we considered no other food source for mussels beside phytoplankton, a strategy that has been successful in previous models , Saraiva et al. 2012

Model domain
Specific information is necessary to accurately describe the model domain and establish a realistic initial point from which to run different simulations. The following information is crucial: • Hydrodynamics -Marine ecosystems are dynamic environments due to winds and tides, and therefore the biogeochemical processes within them are influenced by the circulation of the system. Water exchange coefficients between the different regions of the model domain, far field and river describe the transport processes. A separate physical model is generated based on Navier-Stokes equations, within a grid system, e.g. finite elements. A fully spatial physical model provides a high spatial resolution and, therefore, a more accurate description of the hydrodynamics. Digital coastline and bathymetry as well as sea level data are the minimum inputs to a numerical circulation model. The hydrodynamic model must be independently validated using current direction and velocity, sea level, temperature and/ or salinity. At a coarser level, the estimation of exchange coefficients between connected boxes can be achieved by using temperature time series collected with dataloggers in each box (Dowd 2005) or gradients of natural conservative or semi-conservative tracers (Filgueira et al. 2010).
• River flow -Given the effect of seasonality in river discharge, a detailed time series is required to quantify the contribution of riverine nutrients to the studied area. In addition to gauged stations in the river, GIS tools, such as RiverTools software (rivix. com), can be used to estimate river flow (Peckham 2009).
• Aquaculture activity -Detailed information on aquaculture location, stocking and seeding/harvest schedules is required to accurately allocate biomass throughout the model domain and in time. In addition, information about the population structure is desirable to establish year classes and improve the accuracy of total biomass and its distribution.
• Dissolved nutrients, chlorophyll -Concentrations of these variables are required to define the state of the model domain at the beginning of the simulation.
• Primary production (PP) -PP is the production of the organic compounds derived from photosynthesis and therefore informs phytoplankton dynamics. Time series of PP can be obtained by field sampling or satellite remote sensing (Behrenfeld & Falkowski 1997). The latter suffers from low spatial resolution and poor coverage in nearshore areas. PP could be also modelled, in which case time series of light and turbidity would also be necessary.

Groundtruthing
The correspondence between modelled and observed values must be analyzed to validate model performance. More time series data generally result in more robust groundtruthing. Ideally, components of an ecological model should be groundtruthed in multiple areas of the system. The groundtruthing process has a focus on the following: • Bivalve growth -The use of times series of dry meat content in different locations is the best way to carry out spatial groundtruthing. Bivalves occupy the highest trophic level in these models and their weight integrates the performance of the other submodels. In addition, time series of meat weight provide information that can be used by farmers for management purposes. Bivalve growth time series must be collected by field sampling.
• Dissolved nutrients, chlorophyll -These time series can be used to groundtruth the model when bivalve growth is not available. Although the nutrient and phytoplankton submodels do not represent the highest trophic level, the bivalves interact with these components. Remote sensing tools are not suitable for groundtruthing in small bays where resolution is limited. The alternative in these cases is field sampling.
Generic core for physical−biogeochemical modelling Filgueira et al. (2012) described the steps for creating a generic and flexible core for physical− biogeochemical modelling based on Simile (available upon request). The biogeochemical model can be modified for specific study by means of Simile's user-friendly graphical interface. This generic core ( Fig. 4 in Filgueira et al. 2012) is fully functional for use with multi-box models, where the spatial resolution is defined by large boxes that simulate different parts of the water body and the hydrodynamics are modelled by discrete water exchange between these boxes. However, the authors also described a way to couple the biogeochemical core constructed in Simile with the results of a fully spatial hydrodynamic model (available upon request). Therefore, this gene ric core is able to assimilate hydrodynamic information in 2 different ways to construct either box or fully spatial models, depending on the purpose of the study. This core has been successfully applied in several studies (Grant et al. 2007, Filgueira & Grant 2009, Filgueira et al. 2010 and coupled so far to 3 hydrodynamic modelling packages; AquaDyn, RMA and FVCOM.

Optimization tools
The optimization process is carried out using PEST, an optimization tool that is included in Simile. PEST uses the Gauss-Marquardt-Levenberg algorithm to estimate the value of a parameter or a set of parameters, thereby minimizing the discrepancies between the model results and a dataset chosen by the user. PEST can be used with 2 objectives in mind: (1) tune or estimate parameters in order to calibrate the model, avoiding 'eyeball' calibrations, and (2) optimize management strategies according to a predefined criterion. However, when a dataset is used to tune or estimate parameters, a different one must be used for groundtruthing, otherwise the tuning and the validation are not independent, which would invalidate the groundtruthing. Filgueira et al. (2010) used both capabilities of PEST in a box model to estimate water exchange between different parts of a fjord and the optimum biomass of mussels that could be cultivated according to a carrying-capacity criterion based on chlorophyll depletion. R. Filgueira (unpubl. data) used PEST in a fully spatial model to estimate the best location for mussel culture as well as the optimal mussel density according to a carrying-capacity criterion based on chlorophyll depletion.  Fig. 2A). It is characterized by a mean tidal height of 0.85 m and a maximum tide of about 1.6 m. Depths generally range from intertidal to 19 m, with a single deeper hole (27 m) just inside the inlet (Fig. 3). The centre of the harbour tends to be deepest, with shallower depths in all of the arms. Its tidal prism (the volume of water exchanged on a mean tide) is only 2.37 × 10 7 m 3 , compared to a total estuarine volume (high tide) of 2.87 × 10 8 m 3 , or 8.2% of the total estuarine volume. These characteristics lead to a bulk flushing time of 137 h (~6 d) for St. Ann's Harbour.
Currently, there are 4 operating mussel leases in the harbour ( Fig. 2A), although the 'old leases' were used until recently and are considered in the model. The available register of mussel biomass prepared by the growers since 2003 in the 5 leases provides de -  (Table 2). Two cohorts were considered in the model, which were initialized in order to mimic the situation in the harbour on 1 May as follows: (1) Year 1 mussels with 200 mg dry weight and (2) Year 2 mussels with 800 mg dry weight (R. Stuart pers. obs.). At the beginning of the simulation period the biomass in the estuary was approximately 60 t based on farmers' reports (R. Stuart pers. obs.). The distribution of mussel biomass within the bay was calculated using the final biomass during the harvesting season as described in Table 2. Defining distributions has been challenging, because there have been changes in stocked biomass among the different leases since 2003. For example, years without mussels in a particular lease have not been considered in the calculation of the average biomass in that lease, which could overestimate the total standing stock biomass of the bay. Consequently, the use of this scenario (Table 2) provides a worst-case scenario for analyzing the aquaculture impact on the environment, i.e. the maximum allocated mussel biomass. Available datasets allowed for constructing time series for the period from 1 May until 11 November (195 d), covering the ice-free cultivation period. It was assumed that seeding and harvesting activities exerted a negligible impact on the total biomass allocated in the bay during the modeled period (R. Stuart pers. obs.).
Far field (St. Ann's Harbour) The 8 d frequency time series of 9 km MODIS-Aqua chlorophyll a and monthly 11 µm band nighttime sea-surface temperature (SST) values, averaged within a region located just outside of St. Ann's Harbour defined by the coordinates 46.492° to 46.359°N and −60.243° to −60.409°E, were constructed for the period of 2003 through 2009. These 2 time series were produced with the Giovanni online data system, developed and maintained by the NASA GES DISC. Giovanni is a valuable resource for ecological modelers, as it provides a simple web-based environment to access, explore, analyze and visualize numerous types of Earth science remote sensing data (Acker & Leptoukh 2007).
Although land-based nutrients can be modelled from land use with tools such as N-SPECT (www. csc. noaa. gov/digitalcoast/tools/nspect/index.html), dissolved nutrients derived from the coastal ocean require water sampling or a nutrient climatology. In this case, data from the Atlantic Zone Monitoring Program (www.meds-sdmm.dfo-mpo. gc.ca/ isdm-gdsi/ azmp-pmza/index-eng.html) included a station on their Cabot Strait Line (~60 km from St. Ann's Bay). Data from several years were compiled to generate a time series of nitrate (10 datapoints for the modelled period). In several years of sampling, we have never found detectable ammonia in the bay, and this nitrogen source was not considered further.
There are no available time series of nutrients for the rivers that discharge into St. Ann's Harbour. The land use in the watershed is predominantly subboreal forest on rocky soils. Therefore, it is likely that A 2-dimensional hydrodynamic model constructed with AquaDyn (Hatch Ltd.) is available for the area (J. Grant unpubl. data). The model is based on a spatial grid of finite elements (Fig. 2B), a set of triangular cells in which water volume is tallied at each tidal stage. In the present model, water flux was due to baroclinic circulation only. Turbulence, wind waves and river input were not included, making the model conservative with respect to flushing. Simulations were run for up to 45 d, with time steps from 1200 to 3600 s, the last 29.53 d (lunar month) were considered in the coupling process.
In order to characterize the initial spatial distribution of total nitrogen, chlorophyll and particulate orga nic matter within the harbour, spatial sampling was carried out on 22 July 2010. Fourteen stations were sampled (Fig. 2C), covering rivers that empty into the harbour as well as the boundaries of the model domain. Total nitrogen was measured using a Timberline TL-550 ammonia/nitrate analyzer (Mansell et al. 2000). Two replicates were collected at each sampling point in 25 ml pre-washed polyethylene bottles and frozen (−20°C) until analysis. Water samples for chlorophyll analysis were collected in duplicate. One litre of each replicate was filtered through 25 mm Whatman GF/F filters and kept frozen (−20°C) until chlorophyll extraction. After ex trac tion (90% acetone), chlorophyll was measured fluorometrically (Welschmeyer 1994). Total particular matter (TPM) and particulate organic matter (POM) were measured gravimetrically on pre-ashed (450°C, 4 h), 25 mm Whatman GF/C filters. Two re plicates were collected at each sampling point. One litre of each replicate was filtered, and salts were eliminated by washing with 100 ml of an isotonic solution of ammonium formate (0.5 M). Subsequently, the filters were dried at 110°C for 24 h and weighed to determine the TPM. POM was determined after ashing the filters for 4 h at 450°C.
There were no available direct measurements of PP in the study area. The closest location with available time series of PP was Tracadie Bay (mentioned above). Both bays share the same latitude, which is crucial for determining the daily light cycle and, therefore, PP. Satellite remote sensing was used to establish a correlation between both bays. A monthly time series of the mean net PP was constructed for 2003 to 2009 using satellite imagery. First, monthly averages of global 9 km net primary productivity imagery, based on the vertically generalized production model, were obtained from the ocean productivity website (Behrenfeld & Falkowski 1997, www. science. oregonstate.edu/ocean.productivity/index. php). The imagery was then sub-scened to both areas of interest, Tracadie Bay and St. Ann's Harbour. Finally, the mean value of each monthly subscened image was computed. Both satellite-generated time series followed the same pattern, and a correction function was calculated using satellitegenerated time series to scale the measurements in Tracadie Bay to those in St. Ann's Harbour. A sensitivity test showed that a variation of ±10% in PP causes an average change of + 2.57% / −2.55% in the results of the phytoplankton submodel.

Groundtruthing (St. Ann's Harbour)
The overall correspondence between observed and model values was evaluated by comparing the average observed mussel growth in Year 1 and Year 2 classes with modelled estimations. Two available datasets collected by farmers in all leases of the harbour during 2006 and 2010 (R. Stuart pers. obs.) were used to validate mussel growth estimations.

Physical−biogeochemical model (St. Ann's Harbour)
The results of the hydrodynamic model constructed in AquaDyn were coupled to the biogeochemical one following Filgueira et al. (2012). The biogeochemical model is based on a classical PNZ model (Kremer & Nixon [1978]: phytoplankton [P]−nutrients [N]−zooplankton [Z]), with the addition of mussel (M) and detritus (D) submodels. However, given the absence of zooplankton data and our previous experience that this compartment is of minimal importance (Grant et al. 2008, Filgueira & Grant 2009), zooplankton was not modelled in our case study. The model is characterized in terms of milligrams of C per cubic metre, with the exception of dissolved nutrients, which are ex pressed as milligrams of N per cubic metre. A brief description of the different equations is given in Table 3 with more details in Grant et al. (1993Grant et al. ( , 2007Grant et al. ( , 2008, Dowd (1997Dowd ( , 2005 and Filgueira & Grant (2009). The differential equations are as follows: (1) The model was forced with daily time series, but internally the model was run with a time step of 0.0001 d. Given that some datasets were not col-lected at the same time resolution, linear interpolation was used to calculate daily values.
Optimization procedures (St. Ann's Harbour) The horizontal water exchange between the model domain and the far field was estimated using PEST, performing a calibration between the results obtained in AquaDyn and Simile. Both Simile and AquaDyn were set up in the same way to run a version of the model where a conservative tracer is the only component. Assuming a constant concentration of the conservative tracer at the model boundary (10 U m −3 ) and a homogeneous initial distribution of the tracer inside the harbour at Time 0 (1 U m −3 ), the model was run until equilibrium was reached. PEST minimized the discrepancies between the curves of tracer concentration at different points of the harbour in both approaches by adjusting the water exchange coefficient between the model domain and the far field in a coupled model constructed in Simile.
An optimization procedure was also defined to estimate the carrying capacity, which was defined as the standing stock of mussel biomass that maintained chlorophyll concentration in the harbour within the bounds of natural variation (Filgueira & Grant 2009). The mussel density in the different leases that fulfills this criterion was estimated using PEST by minimizing the discrepancies between the chlorophyll time series observed in an aquaculture scenario and the time series observed in a theoretical background  Table 3. General model terms and sources for detailed explanations scenario without aquaculture. These background chlorophyll time series were adjusted to represent the lower levels of chlorophyll that include the range of natural variation observed in the harbour. In this way, the mussel biomass at the calculated carryingcapacity level may have reduced chlorophyll concentration in the harbour, but never beyond the limit of its natural variation. The natural variation of chlorophyll in St Ann's Harbour was obtained by analyzing a 7 yr monthly chlorophyll time series observed in the far field. The highest and lowest chlorophyll values in each season were removed from the analysis to eliminate potential outliers. This approach is a conservative way of analyzing natural variation by reducing the dispersion of the data while maintaining average values. Analysis of the time series yielded a minimum seasonal Pearson's coefficient of variation of 32.3% (Table 4). This value is considered to be the natural variation limit for the carrying-capacity aquaculture scenario optimized by PEST.

Conservative tracer simulations
The tracer simulations in both AquaDyn and Simile eventually resulted in a tracer concentration equilibrium between the bay and the far field, but prior to that time, the distribution of tracer could be interpreted as a proxy for flushing. The tracer concentrations after 25 d of simulation of both AquaDyn and Simile models indicated that there was a strong gradient from the mouth to the inner harbour (Fig. 4). The penetration of the tracer was slightly asymmetrical, and the average conditions used in Simile indicated that there was reduced exchange on the southern side of the harbour compared to the northern side. These results demonstrated that Lease 1188 had the highest flushing, followed by 1187, the old leases, 1189, and 1186, respectively. The optimization process carried out to determine the instantaneous water exchange between the model and the far field yielded a value of 17.25 m 3 s −1 . A sensitivity test carried out on this estimated water exchange parameter of ±10% resulted in a change of + 3.70% / −3.93% in the average tracer concentration in the entire harbour after 25 d.

Groundtruthing
Mussel growth simulated by the model was within the range of variation observed in the bay in both year classes (Fig. 5). Mussel dry weight estimates for the Year 1 class were within the confidence limits of the observations with the exception of 1 datapoint on 13 July, which was higher than the modelled values as well as the other field observations. Modelled values for Year 2 mussels were also within the confidence limits of the observations, with the exception of the sampling carried out on 30 September. The observed mussel weight on 30 September indicated negative growth during September which might have been caused by spawning. The simplified model used in this study did not include spawning, and, consequently, the model could not explain the growth pattern observed in September.

Implications of stocking biomass for mussel growth and chlorophyll depletion
Scenario building was used to determine the effect of variations in stocking biomass on mussel growth and culture yield. The evolution of average individual growth within the bay is represented in Fig. 6 for 6 different scenarios, with initial stocking biomass ranging from 15 to 90 t. At the start of these simulations all scenarios showed the same performance. However, after 50 to 75 d (19 June to 14 July), the individual growth of both the first and second year classes became density dependent, i.e. reduced growth for higher stocking densities. For example, in Year 2 mussels there was a growth penalty in individual weight of up to 22.5% for the 90 t scenario compared to the 15 t scenario. This effect was even greater for Year 1 mussels, reaching a penalty of 28.1% for this same stocking comparison. This im pact on individual growth  Table 4. Averaged seasonal chlorophyll concentration and natural variation analysis. CV: coefficient of variation reflected the limited capacity of the environment to provide food and support the growth of an entire population. The evolution of averaged chlorophyll depletion within the harbour is represented in Fig. 7 for the same 6 scenarios. The percentage of chlorophyll de pletion increased through time due to the corresponding increase in mussel biomass. Depletion was compared to the minimum seasonal coefficient of variation in chlorophyll in the far field, 32.3%. Chlorophyll reduction by suspension-feeders below this threshold (32.3%) was within the natural variation of the system. The different scenarios were grouped into 2 classes: low initial stocking scenarios of 15, 30 and 45 t, which never caused a depletion in chlorophyll beyond the natural variation threshold, and high initial stocking scenarios of 60, 75 and 90 t, which passed this threshold on 25 September, 21 August and 24 July, respectively. The current scenario in St. Ann's was estimated to have a stocking biomass at the beginning of the simulation of 60 t. Therefore, depletion would be below the natural variation threshold for about 75% of the studied period.

Chlorophyll depletion at constant mussel standing stock biomass
Changes in stocked biomass via growth and variation in boundary conditions may combine to cause extreme values in chlorophyll depletion. In order to estimate a steady state chlorophyll depletion in the harbour, new simulations were run with constant mussel biomass in the whole harbour. This approach . Time series of chlorophyll depletion for 6 different scenarios with initial stocking biomass of mussels ranging from 15 to 90 t. The horizontal grey line represents the sustainability threshold (32.3%) based on natural variation analysis assumes that mussel biomass interacts with the ecosystem model as a forcing function rather than a response variable (Dowd 2005). The mussel compartment was fully functional in the model, but its biomass remained constant through time. Following this approach, maps of averaged chlorophyll depletion through time were generated for several constant standing stock scenarios. For a constant standing stock of 60 t, average chlorophyll at bay-scale was reduced by ~24.9% (Fig. 8A). Doubling of this density yielded bay-scale chlorophyll depletion levels in the range of ~41.2% (Fig. 8B). For the very high standing stock density (180 t), bay-scale chlorophyll was severely reduced to ~51.7% (Fig. 8C). These results can be summarized in a curve of average depletion versus cultured biomass (Fig. 9). With increasing crop size, chlorophyll depletion reached an asymptote. It should be noted that the asymptote was < 75% depletion since all of the water and its chlorophyll was not accessible to the farm sites. In some senses this was a worst-case scenario, since mussel biomass was continuously high rather than building through growth. The current production scenario was estimated to have a stocking biomass of 60 t at the beginning of the simulation, which produced an averaged maximum depletion of 24.9%, below the natural variation threshold. The worst-case scenario according to actual aquaculture activity (Table 2) resulted in an averaged standing stock biomass at harvest of 191.3 t, and an averaged chlorophyll depletion of 53.3%, beyond the natural variation threshold. Following the current distribution of biomass within the bay (Table 2) and according to the results in Fig. 9, the natural variation threshold would be reached with a constant standing stock biomass of 82.2 t.

Carrying-capacity estimation: sustainable aquaculture scenario
An optimization procedure carried out using PEST was performed to estimate the standing stock biomass that fulfills the carrying-capacity criterion (depletion within natural chlorophyll variation). The averaged standing stock biomass predicted by PEST was 84.5 t and the spatial distribution was 1.35, 1.71, 2.29, 0.00 and 2.89 g m −3 in the old leases, and Leases 1186, 1187, 1188 and 1189, respectively. This distribution implied a significant change from current practices because it recommended that culti vation of mussels in Lease 1188 be avoided. The spatial distribution of averaged chlorophyll depletion through time is presented in Fig. 10, and the averaged value for the whole harbour was 32.3%. The optimum standing stock biomass predicted by PEST was close to the calculation derived from the current distribution of biomass within the bay (Table 2) and the results in Fig. 9, 82.2 t.

Actual scenario with population dynamics
According to the current scenario, 60 t at the beginning of the simulation following a proportional allo-cation according to Table 2, the cultured biomass reached the natural variation threshold on 25 September (Fig. 11). This indicated that chlorophyll depletion in the harbour was higher than the natural variation for ~25% of the year. This may have been caused by the absence of population dynamics in the model design. There were no local data to initiate a population dynamic model, so a simple component was added to the main model in order to explore the consequences of mortality. Following Mallet & Carver (1989), a mortality rate of 20% yr −1 was used for both year classes. The 60 t scenario using these mortality rates resulted in a maximum chlorophyll depletion lower than the natural variation threshold (Fig. 11), reducing the depletion by ~7% by the end of the simulated period compared to the scenario without mortality. Thus, mortality may be an important factor for evaluating the potential ecosystem implications of mussel aquaculture.

Sustainable mussel culture
Maintaining the percentage of phytoplankton depletion within the bounds of natural variation is a straightforward way to manage sustainability from an ecosystem standpoint. Given that phytoplankton are the base of marine food webs and that filterfeeder populations are important for their control (Cloern 1982, Dame & Prins 1997, the stress that bivalve aquaculture exerts on the ecosystem has been quantified in terms of seston depletion (Campbell & Newell 1998, Pouvreau et al. 2000. Average chlorophyll depletion for the optimal carrying-capacity scenario according to optimization by model-independent parameter estimation (PEST, www. pesthomepage. org). Lease areas outlined and lease numbers shown. No mussels were allocated in lease 1188 (grey) in PEST scenario Fig. 11. Time series of chlorophyll depletion in the 60 t scenario with and without mussel mortality. The horizontal grey line represents the sustainability threshold (32.3%) based on natural variation analysis , Ferreira et al. 2007, Duarte et al. 2008). We propose the use of ecological resilience as a framework for quantifying the limits beyond which ecosystem health is compromised (Fig. 1). We assume that the natural variation of ecosystem variables is below tipping points beyond which the stability of the system is endangered. Consequently, managing within natural variation thresholds provides precautionary limits that are capable of supporting unexpected disturbance without reaching the tipping points and, therefore, of guaranteeing ecosystem sustainability.
According to the simulations that included mortality rates, the current aquaculture level in St. Ann's Harbour resulted in maximum chlorophyll depletion of 29.6%, which was within the threshold of natural variation, 32.3%. The optimization procedure applied using PEST to calculate standing stock biomass that fulfills the carrying-capacity criterion provided a different distribution of biomass in the leases compared to the actual one. Re-distribution of the biomass according to the optimized solution allowed an overall increase of 0.3 t. The change in spatial biomass allocation suggested by PEST did not substantially change chlorophyll depletion compared to the existing farm arrangement. Therefore, the actual mussel distribution among the leases seems acceptable from a sustainability perspective. However, given the limitations in model parameterization and validation derived from the use of data with high uncertainty, such as the nutrient time series (ocean and river) or scarce groundtruthing data, the results generated for St. Ann's Harbour cannot be considered a final and static outcome. On the contrary, the results should be understood as the best objective scientific assessment that is possible to create with the available data. Further experiments must be performed to reduce uncertainty in forcing time series, as well as to carry out a full and rigorous validation of all the components of the model. St. Ann's Harbour was chosen as a very challenging example of a data-poor environment and, consequently, was an ideal testing ground for the techniques described in this study. The case study is considered as proof of concept; however, calculations should be understood as initial assessments of carrying capacity rather than conclusions for decision-making.

Collecting datasets and uncertainty analysis
Modelling exercises demand an extensive amount of data. In this work we emphasized 3 aspects re -lated to data collection: collaboration with industry, generation of time series using remote-sensing techniques and collation of existing sources. Collaboration with stakeholders is vital for marine spatial planning, for instance, in the case of farms providing the information needed to resolve biomass distribution within the bay. In addition, protocols for tracking biomass allocation, seeding and production can facilitate the use of datasets in future research projects. For example, the 7 yr time series used in this study to estimate the mussel biomass distribution throughout an entire bay is an exceptional dataset for use in aquaculture-related modelling exercises. The absence of these datasets would increase the cost of further research, as well as the uncertainty of the results.
Satellite remote sensing data are particularly advantageous for modelling exercises, since these data allow consistent synoptic observations of a study area for periods ranging from years to decades. In addition, the data are typically free and easily accessible, and many tools exist to facilitate analysis and visualization (e.g. Giovanni). There are limitations, however, to the use of satellite remote sensing data, particularly regarding the resolution in coastal regions. For example, MODIS Aqua spatial resolution of chlorophyll is at best 500 m per pixel and 250 m per pixel for turbidity, requiring processing of level 0 data (Hu et al. 2004, Chen et al. 2007). In addition, although satellite remote sensing is widely considered to be successful in the open ocean, particularly regarding biological parameters, in coastal waters that are typically more optically complex there has been variable success. This is due to the fact that the algorithms used to retrieve water constituents from satellite re mote sensing data, as well as essential atmospheric corrections, are not necessarily applicable to coastal waters (Moses et al. 2009).
The use of historical datasets complemented with sensitivity tests to evaluate the impact of source variability on the model's outputs could be a solution in data-poor environments. Another possible tool for application in data-poor environments is the use of 'extreme scenarios'. For example, given the absence of nutrient time series in the main river that empties into the harbour, values collected in a nearby location, but one with intensive agricultural activity, can be used instead. This model configuration provides an extreme situation with regards to the land-use pattern in the watershed of St. Ann's Harbour. However, it allows the analysis of an 'extreme scenario' and the quantification of the maximum uncertainty of the modelling exercise.

Ecological modelling for ecosystem-based management
Ecological modelling has been successfully applied to aquaculture sites to evaluate the potential of new cultivation areas (e.g. Brigolin et al. 2006, Filgueira et al. 2010, optimize the profitability of existing farms (e.g. Héral 1993, Ferreira et al. 2009) and manage ecosystem-level impacts (e.g. Grangeré et al. 2010, Byron et al. 2011). These models vary in complexity from simple ratio-based budgets (e.g. Dame & Prins 1997) and box models (Pastres et al. 2001) to sophisticated fully spatial hydrodynamic−biogeochemical coupled models (e.g. Ferreira et al. 2008). The approach described in the present paper combines existing models, historical time series and remote sensing datasets into a fully spatial model to provide a management strategy for optimizing aquaculture production from a sustainable ecosystem standpoint, according to FAO recommendations (Soto et al. 2008). Similarly, Silva et al. (2011) have carried out a study to evaluate the best site for shellfish aquaculture in Valdivia (Chile), which combines GIS-based models with a farm-scale carrying-capacity model (Ferreira et al. 2007). As in our work, Silva et al. (2011) developed a modelling exercise for a study site with limited data availability, simplifying the model according to available datasets. Despite the simplification, Silva's approach is focused on critical dynamics of the domain area, consequently providing valuable results for decision-makers in the practical application of ecosystem-based aquaculture. A simplification in the model design could optimize the trade-offs between the effort needed to carry out the modelling exercise and the uncertainty of model outcomes (Nobre et al. 2010).
Anticipating future conditions is an important goal of ecosystem-based management and marine spatial planning (Polasky et al. 2011). Ecosystem modelling is the perfect tool to anticipate future conditions, but implementation of ecosystem models has many data requirements, which are both time consuming and costly. Our study reviews the minimum datasets that are needed to construct a model to evaluate the impacts of bivalve aquaculture in data-poor environments, as well as alternative techniques for collecting these datasets. Tools such as ecoinformatics, remote sensing, scenario analysis and optimization can help fill these data gaps. However, the use of these techniques increases the uncertainty of model results, and, consequently, straightforward application to decision-making is not trivial. Model outcomes in data-poor environments should be understood as the best objective scientific assessment that is possible to provide with the available data, which is crucial in the field of applied sciences (Polasky et al. 2011). This is the case for St. Ann's Harbour, which has been chosen as a testing ground for the techniques described in this study. Therefore, the results of the St. Ann case study must be considered in the context of high uncertainty as the best scientific assessment of carrying capacity rather than a final outcome for decision-making. In addition, the results could be used to provide feedback for planning processes via adaptive management (Polasky et al. 2011, Halpern et al. 2012, which relies on the understanding of marine spatial planning as a dynamic process rather than a static outcome.