NUMERICAL MODELING AS SUPPORTING TOOL FOR AQUACULTURE OF OYSTERS IN A SUBTROPICAL ESTUARINE ECOSYSTEM

The Estuarine-Lagoon Complex of Cananeia-Iguape (ELCCI) is a large estuarine system dominated by mangroves, located on the southeast coast of Brazil (25° S, 48° W). Oyster farming is an important economic activity in Cananeia due to the natural abundance of oyster seeds originally found in its mangrove fields. The aim of this study was to describe the hydrodynamic aspects of the mangrove and oyster farming in the ELCCI by analyzing environmental scenarios using numerical simulations. Simulations allowed for the description of the distance of penetration of the salt wedge into the system as well as the direction and magnitude of the tidal currents. Higher values of river discharge corresponded to lower distances of penetration of the salt wedge, a feature which shows that the freshwater discharge and geomorphology together determine the variation in salinity within the estuary. Oysters are cultivated in regions of higher salinity and relatively stable conditions, with low speeds of tidal currents and little vertical variation in salinity. Under average river discharge conditions (450 m3.s-1), there is a confluence of tidal waves in Mar Pequeno, where velocity is lower. This change may lead to a higher sedimentation at restricted areas and thus to a higher deposition of organic matter derived from oyster cultivations during flood tides.


INTRODUCTION
Estuaries are highly dynamic environments and are subject to constant environmental changes -both natural and anthropogenic changes. The changes in these systems are forced by factors such as solar radiation, tides, winds, precipitation, and river discharge (Moura and Nunes, 2016). Tides define the barotropic component of the pressure force that drives circulation, which induces most of the horizontal movements. Differences 2/10 in density promote the baroclinic component that pushes the denser water along the bottom, up the estuary defining the salt wedge. Altogether, tides, river runoff, stratification, and stirring processes are the factors that define the estuarine circulation patterns (Miranda et al., 2002), as well as saltwater intrusion and the salt wedge positioning. Salinity is an important environmental factor for determining estuarine oyster distribution.
The Cananeia estuary is considered one of the most important wetlands on the Brazilian coast in terms of biodiversity and natural productivity, being recognized both nationally and internationally as the third most productive ecosystem of the South Atlantic, due to its well-preserved characteristics (UNESCO, 2005). It possesses favorable conditions for the creation of natural banks as well as for the implantation of farms for the fattening and harvesting of the oyster Crassostrea sp.. This region is responsible for the largest oyster production in the Brazilian state of São Paulo (Pereira et al., 2003).
According to Miyao and Harari (1989) salinity is the key parameter in determining patterns of estuarine circulation. Salinity has a positive correlation with tidal range, and it allows for the calculation of exchange rates between estuary and adjacent oceanic regions when temperature gradients are not large enough. Salinity should be considered in the malacoculture, since it presents daily and seasonal variations in the estuaries, being influenced by the tidal regime and the rainy season (Vilanova and Chaves, 1988;Funo et al., 2015). Moreover, it is an important environmental factor that determines the distribution of bivalves in estuarine and marine environments (Fürsich, 1993).
Studies have shown that the environmental conditions strongly impact the survival of bivalve mollusks during some phases of their life cycle (Cáceres-Puig et al., 2007;Dickinson et al., 2012;Funo et al., 2015). These conditions are related to factors such as temperature, salinity, pH, presence of microalgae and composition of particulate matter in suspension (Guzmán-Agüero et al., 2013).
Numerical modeling of hydrodynamics in estuaries can be an innovative tool to support farming of oyster and even other aquatic organisms in estuarine ecosystems. Especially in the studied case of the ELCCI, the region is very suited to attend the demands of the fish producing market since there are plenty of conditions supported by the local infrastructure and community tradition on fishery and food processing.
The objective of this study was to elucidate the best location choices for oyster farming by applying the Delft-3D suite of models to simulate the hydrodynamics of the Cananeia-Iguape estuarine system. Modeling specific environmental scenarios allowed us to describe the oyster farmings susceptibility to the hydrodynamic aspects of the Cananeia-Iguape estuary by defining the salt wedge positioning in the estuary from the salinity distribution in the system. The climatologic and oceanographic conditions defined the most significant scenarios on the estuarine dynamics. Finally, modeling results revealed the general circulation pattern in the estuary that tends to favor the best locations for mangrove setting for oyster farms.

Study area
Located on the southern coast of São Paulo state (25° S, 48° W), the ELCCI is a region dominated by mangroves and extensive areas of rainforest. As seen in Figure 1, the ELCCI is 75 km long and consists of four large islands and six main channels: Baia Trapande, between Ilha do Cardoso and Ilha de Cananeia; Mar de Itapitangui and Mar de Cubatão, between Ilha de Cananeia and the mainland; the Mar de Cananeia, between Ilha Comprida and Ilha de Cananeia; Mar Pequeno, between Ilha Comprida and the mainland; and Valo Grande channel, a man-made passage between Ilha de Iguape and the mainland. The system communicates with the Atlantic Ocean at two entrances: the Barra de Cananeia (South), and the Barra de Icapara (North), which receives freshwater and sediment from Rio Ribeira de Iguape. The southernmost portions of the complex are also influenced by small rivers, which form a drainage basin of 1,339 km 2 (Mishima et al., 1985). The Valo Grande channel, that connects the Ribeira de Iguape River to Mar Pequeno, is responsible for diverting 70% of the fresh water from the river into the estuary.
Tides in the region are predominantly semi-diurnal, with mean amplitude of 83 cm during spring tide and 13 cm during neap tides (Miyao and Harari, 1989). The main active components of the tide in the study area are O 1 , M 2 , K 1 and S 2 . They are the more energetic tidal components and therefore, the most important components for defining the local tides (Harari and Camargo, 1994). Mangrove oyster farmings are predominantly found nearby Mar de Cubatão, Baia Trapande and Mar de Itapitangui where the environment is more stable and the influence of small rivers allows for the development of large mangroves.

Numerical modeling
The Delf3D Model solves the three-dimensional Navier-Stokes equations for an incompressible fluid, with the non-hydrostatic and Boussinesq approximations for a non-inertial system of reference. Details on the formulation and discretization of momentum equations can be found in Deltares (2013). The hydrodynamic module employs a horizontal, curvilinear and orthogonal mesh, capable of best fitting the geographical representation of the ELCCI region into the computational domain. Vertical discretization uses a sigma coordinate transformation, which ensures that depth variations maintain a proper resolution into both the superficial and bottom boundary layers. Indeed, it is a four-dimensional model, because it allows the simulation to run within the time domain, an option which best evaluates the estuary dynamics of tides, river intakes, mixing processes, and buoyancy effects. The system of equations includes the diffusion of heat and salt, turbulence, continuity and a state equation for density.
The ELCCI computational domain, created on Delft-3D-RGFGRID module, consists of five different grids, each of which resolves the dynamics of coastal and estuarine regions with proper space and temporal resolution, respectively terms of magnitude O(10 1 to 10 2 s) and O(10 1 to 10 3 m), as in Figure 2 and Table 1. The bathymetry of the study area was acquired by digitizing the Brazilian Navy's nautical charts N. 1702 and 1703 (DHN MM 1952, MM-DHN 1985. The Delft-3D-QUICKIN module included the Delaunay triangular interpolation method for allowing a proper numerical interpolation of the bathymetry ( Figure 1).
The time running of simulations considered the period of May 7, 2012 to May 29, 2012, with a time step of 30 seconds, for all five domains: valog, iguape, sea1, sea2 and coast. Details of domains are summarized in Table 2 and modeling conditions in Table 1. The model started its calculation from a resting state, so the simulations were run for five days to achieve hydrodynamic balance and computational stability of the system. The model was then adjusted to present hourly results, during a period that corresponded to a complete fortnightly spring-neap tidal cycle, in order to retain most of hydrodynamical features related to tides. We neglected the influence of winds based on the assumptions made by Geyer (1997), in which the influence of winds on the dynamics of shallow and mouth constricted estuaries can only be considered during strong wind events. The spin-up time for model initialization was 48 h, suffice to assure that transient startup noise was dissipated.

Model forces and boundary conditions
We defined three river discharge values for the Rio Ribeira de Iguape, which were extracted from rating curve tables (São Paulo, 2008). This database is constantly updated since its beginning in the late 1950's. The time of response of the estuarine system relates the river discharges with their probability of occurrence over time. Discharge values were used in the following simulations:  • Q 95 : Minimum discharge, which occurs 95% of the time and corresponds to the dry season of the river; • Q 5 : Maximum discharge, which occurs 5% of the time and corresponds to the period of river flooding; • Q M : Long-Term Average of the discharge, which represents the yearly average discharge.
Salinity, temperature and currents data used in the modeling were acquired from the literature (Bernardes and Miranda, 2001). For temperature and salinity of the coastal domain (coast), we used in situ measurements, which were collected on October 7, 2011. Table 2 shows the boundary conditions for each grid used in the implementation of the model.
There are two connections of the estuarine system with the open ocean. Thus, we defined the sea level boundary conditions as Barra de Cananeia (South) and Barra de Icapara (North) from where we used data from two gauges (Table 3): Cananeia (FEMAR, 1956) and Iguape (FEMAR, 1982). Only the major diurnal and semidiurnal components were retained, in order to represent the fortnightly modulation of tides.

Simulations
Three simulations were run throughout a whole fortnightly interval to define the spring and neap periods. Each simulation considered the hydrological conditions as previously defined by river inflows.

Simulation 1
Using Q 95 as defined, the results intended to represent the estuarine condition during drought events. Four scenarios were selected within the whole fortnightly period: spring high tide (Q95_SHT); spring low tide (Q95_SLT); neap high tide (Q95_NHT) and neap low tide (Q95_NLT).

Simulation 2
Using Q 5 as defined, the results intended to represent the estuarine condition during the peak of stormy events. Four scenarios were selected within the whole fortnightly period: spring high tide (Q5_SHT); spring low tide (Q5_SLT); neap high tide (Q5_NHT) and neap low tide (Q5_NLT).

Simulation 3
Using Q M as defined, the results intended to represent the estuary during its hydrologic mean condition. As before, the scenarios were selected within the whole fortnightly period: spring high tide (QM_SHT); spring low tide (QM_SLT); neap high tide (QM_NHT) and neap low tide (QM_NLT).

RESULTS
In order to demonstrate the influence of the estuarine hydrodynamics in mangrove oyster farming, we considered the results of salinity and currents in the sea1 domain (from Barra de Cananeia to Barra de Icapara) because it comprises most of Mar Pequeno, throughout the West channel. The sea2 domain (East channel) is equally relevant, following roughly the same dynamic aspects felt by its counterpart, the West channel. For the sake of better comprehension, only relevant simulations were taken, and salinity was described as oligohalyne (0.5 < S < 5.0); mesohaline (5.0 < S < 18.0) or polyhaline (18.0 < S < 30.0).
All simulations showed that waters entering Barra de Cananeia divert in Mar de Trapande, flowing through the West and East channels, and merging again at the single and long path at Mar Pequeno. The Ribeira de Iguape diverts at Barra de Icapara when part of its discharge flows to the coastal region, and the other path flows Westward through Mar Pequeno. Therefore, there is a tendency of tombolo formation somewhere in Mar Pequeno, between Barra de Cananeia and Barra de Iguape. When the water level is high at the bars due to flooding, it is low at Mar Pequeno, and during low tide, flow is reversed.
Mesohaline condition was observed in Mar Pequeno (S in 5-15 range), while polyhaline conditions occurred in Baia Trapande (S > 20). On the banks of Baia Trapande, where the trays of oysters are located, the current velocity is always less than 0.2 m.s -1 , and currents flow unidirectionally during both flood and ebb tides. Salinity had a maximum range of five units between high and low spring tides. From Baia de Trapandé to Mar Pequeno, throughout West Channel, polyhaline waters extend for 35 km with low variation in salinity for both spring and neap tides. During neap tide, this area was partially mixed, while during spring tide, complete vertical mixing dominated. Figure 3 shows the direction and velocity of currents in the estuary under conditions of minimum river discharge, regarding this as a relevant scenario. The merging of tidal waves occurs at Mar Pequeno, where the velocity of the currents varies between 0 and 0.2 m.s -1 . The salinity vertical profiles were taken along the sea-1 domain, from Barra de Cananeia (km 0) to Barra de Iguape and Barra de Icapara (90+ km), with two salt wedges dynamically positioned from the bars. Polyhaline waters entering Barra de Cananeia flows up estuary, so mesohaline waters penetrade more than 60 km during spring tide. The salt wedge from Barra de Iguape and Icapara, however, extends for 10 km during high tide (Figure 4). Simulation 2 (Q 5 ) Figure 5 shows the direction and velocity of high tide currents in the estuary under conditions of maximum river discharge during the flood tide, when water entering in Barra de Icapara flow through Mar Pequeno, in the direction of Cananeia, partially piling up the river outflow. During the ebb tide, the flow continues toward Cananeia Bar, now increasing the river contribution flow. The speed of the current in Cananeia varies between 0.05 and 0.3 m.s -1 , and tidal waves merge at Mar Pequeno (km 30-40).

Simulation 1 (Q 95 )
The salinity wedge entering Barra de Cananeia extends for 35 km, while the water entering Barra de Icapara reaches only 5 km during high tides (distance 90 km, in Figure 6). Freshwater flows through Valo Grande and, in doing so, occupies Mar Pequeno and reaches position 40 km ( Figure 6). Figure 7 shows the direction and speed of high tide currents in the estuary under conditions of average river discharge. Again, tidal waves merge in Mar Pequeno, where weak currents vary between 0 and 0.1 m.s -1 . During quadrature, salinity extends for 60 km from Barra de Cananeia and, during spring tides it extends for almost 70 km (Figure 8).

Model validation
The assessment of modeling projects usually relies on a skill parameterization method, as that proposed by Willmott (1981), in which it evaluates a statistical ranking that measures how much the model results adhere with time series measurements of salinity, temperature, currents, and sea level. Due to the lack     7/10 of consistently oceanographic time series measurements in this project, this evaluation relied on data reported in literature and measurements of temperature and salinity from some CTD campaigns. Hence, we applied a simpler model assessment based on the Pearson correlation coefficient of measurements taken during combined events of high and low waters; syzygy and quadrature and hydrological regimes. The results were C s = 0.64 and C t = 0.73 for salinity and temperature coefficients; C u,v = [0.82, 0.66] for cross and along estuary velocities, and C sl = 0.87 for sea level.
On the model validation, the entire agreement between observed and modeled results will lead to one, while complete disagreement will result in zero. The obtained Results of 0.64 and 0.73 represent a fair agreement between modeled and observed results for the thermohaline variables, while velocity and sea level displacements were the best fit.

Surface circulation
The results point out to a greater marine influence in the southernmost portions of the system, where there are intense currents (Figure 3). This difference should be associated with the shape of Barra de Cananeia, which is perpendicular to the shoreline, thus facilitating the entry of water during the flood tide (Tessler and Souza, 1998;Barbieri et al., 2014). Furthermore, the tidal wave propagates northeastward along the coast (Picarelli et al., 2002), so the tidal currents flow into the system through Barra de Cananeia and not Barra de Icapara.
The longest distance traveled by salt water in the sea1 domain may be associated with the position of Trapande Bay on the coast. According to Paolo and Mahiques (2008), salt water enters preferentially through this bay, and it flows off more rapidly through the Cananeia Sea. Timing delays observed between high tide and low tide at Barra de Cananeia and Mar Pequeno may be partially related to the estuary shape and its morphology, lowering current speeds by energy dissipation through the bottom and lateral stresses. That feature makes the water to take nearly six hours to reach the center of the system. Miranda et al. (1995) explain that these gaps as a result of a combination of progressive and opposite propagating waves, generating a stationary wave in this estuary.
Overall, current velocities (from 0 to 1.0 m.s -1 ) are consistent with previously published work by Miranda et al. (1995) and by Tessler and Souza (1998), who reported that speeds ranged between 0.03 and 1.2 m.s -1 . Miranda and Castro-Filho (1996) measured current speed along the Cananeia Sea about 0.026 m.s -1 during ebb tides. However, the present study could not verify significant differences between ebb and flow tides.
The scenario presented in Simulation 3 (QM = 450 m 3 .s -1 ) represents what most frequently occurs in the Estuarine-Lagoon Complex of Cananeia-Iguape, by consideration of a long-term dataset of hydrologic data. The maps of modeled tidal currents corroborate the results described by Tessler and Souza (1998) and Picarelli et al. (2002).
Results of Simulation 2 (Q5 = 1500 m 3 .s -1 ) show a major influence of freshwater in the estuarine system when compared to the results of other simulations. Maps of current speed obtained in this simulation may probably be correlated with periods of high rainfall in the system, which occurs from January to March (Bergamo, 2000). During that unusual discharge conditions, freshwater from the Ribeira de Iguape River can reach the region where oyster cultivation sites are located.

Salinity
In ELCCI, salinity varies with tide, river discharge and location in the estuary, as shown in the present study. Higher flow values from the Ribeira de Iguape River corresponded to a lower distance of penetration of the salt wedge, a finding that shows that freshwater discharge and geomorphology determine the distribution of salinity in the estuary.
Maps of salinity by Mishima et al. (1985) agree with the results obtained in Simulation 3 (QM = 450 m 3 .s -1 ), in which the freshwater reaches the connection of Mar Pequeno with Cubatão Sea. Salinity levels in the Cananeia Sea also resemble other studies, in which salt water was restricted to Barra de Cananeia, without variation between high tide and low tide. The closer the water is to the Barra de Icapara or Barra de Cananeia, the higher the salinity is, due to the direct influence of the ocean. As it moves away from the river exits, the transport of salt will occur through turbulent diffusion (Miranda et al., 1995). Gavés et al. (2015) observed that the oysters submitted to the salinity of 5 remained with the valves closed in the first two days and only opened them after the third day, coinciding with the increase of mortality of the oysters submitted to this treatment. In the other evaluated treatments, the oysters remained filtering with the valves open and feces deposited in the bottom of the culture containers were observed daily.
Nearby rivers from the Cananeia basin have relatively lower values of salinity but never reach zero. This phenomenon likely occurs because the influence of the sea bringing salty waters into the estuary is more relevant than the freshwater discharge (Coimbra et al., 2007). This factor may explain why the locals carry out the cultivation of oysters in the system once the environment is more hydrodynamically stable, and where salinity promotes the growth of bivalves. This fact may also be associated with the formation of some pools of high salinity at Baia Trapande since these pools promote the growth of oysters and provide places for larvae settlement (Akaboshi and Pereira, 1981). According to La Peyre et al. (2013) salinity less than 5 occurs during periods of lower temperatures (<25 °C), since low salinity and high temperature favors a lethal combination for oysters.
The waters of the Cananeia region may be classified as salty following the CONAMA resolution nº 357 (Brasil, 2005) as they present salinity from 0.5 to 30. This resolution establishes as the norm for the cultivation of bivalve mollusks for human consumption, a pH of between 6.5 and 8.5, and dissolved oxygen (DO) above 5.0 mg.L -1 . The environmental data registered in this study presented values within the parameters defined by the resolution quoted for Class-1 salty water, for the attendance of environmental legislation on conservation and aquaculture.
According to Ferreira et al. (2007), mangrove systems are generally composed of sediments of the silt and clay fraction, due to the low hydrodynamics predominant in these places; however, the water flow velocity is variable in mangrove environments and may cause variations in grain size and distribution. In addition to that, the predominance of fine sediments also reflects the local hydrodynamics. The deposition of fine sediments is not possible in areas with strong currents, for this reason, oysters are hardly found in areas with strong currents.
The environmental impact of aquaculture in an estuarine system can be evaluated by the amount of organic matter that is suspended and deposited on the floors of the rivers and bodies of water (Barbieri et al., 2014). Thus, according to the direction and speed of currents, and according to the patterns of salinity, this deposition will occur at different locations in the estuary and different rates. Furthermore, all sedimentary features of the estuary are directly associated with the dominant currents (Tessler and Furtado, 1983).
Although in this study we did not modeled sediment transport, the slow speeds in Mar Pequeno for all simulations suggest that high rates of sedimentation may occur at that location. Tessler et al. (1990) found high rates of deposition in Mar Pequeno, which corroborates our results on hydrodynamic aspects. According to Barcellos et al. (2009), sediments from the southernmost portions of the system have lower rates of deposition. In the Valo Grande canal, higher sedimentation rates occur due to rapidly decreasing of depth at Mar Pequeno (Mesquita and Harari, 2003;Saito et al., 2006). This change may be the cause of the higher concentration of organic matter at this location. Barcellos et al. (2009) found that the organic matter in the system is well mixed in the estuary, except in Valo Grande canal. Where the tidal waves merge, it is expected that sedimentation rate will be higher, and deposition of organic matter must also occur. Thus, part of the organic matter from the cultivation of oysters must be deposited at this location and must be carried by flood tides and imprisoned by the current from the Cananeia Sea.
Bivalve are filtering organisms which feed on the particles and microalgae which are found in the water and which accumulate in their tissues, including large quantities of organic and inorganic substances, as well as the microorganisms present in the environment, thus acting as bioindicators of unhealthy quality of the water (Pereira et al., 2003;Barros and Barbieri, 2012). It was for this purpose that we studied the sedimentation dynamics because this factor can determine the distribution of oysters in the ELCCI.

CONCLUSIONS
Simulations allowed for the description of the distance of penetration of the salt wedge into the system as well as the direction and magnitude of the tidal currents. Higher values of river discharge corresponded to lower distances of penetration of the salt wedge, a feature which shows that the freshwater discharge and geomorphology together determine the variation in salinity within the estuary. Oysters are cultivated in regions of higher salinity and relatively stable conditions, with low speeds of tidal currents and little vertical variation in salinity. Using maps of currents and salinity distribution in the region, farmer, stockholders, and local fishing communities can benefit from the right choices on the best location for new sites for oyster farming. Also, from climate forecasts, it is possible to know beforehand how the seasonal and extreme events can interfere in the estuarine regime. Model simulations, therefore, can be used to prevent loss of productivity and the occurrence of endangering events that affect aquatic farming activities.