Water quality impact from shrimp farming effluents in a tropical estuary

Coastal ecosystems, especially estuaries, are submitted to environmental impacts caused by excess nutrients from natural processes and anthropogenic activities (Barcellos et al. 2019) in contributing watersheds. In South America, estuaries face uncontrolled occupation on their banks by urban and industrial centers, as well as increasing agriculture and aquaculture activities (Ferreira et al. 2011; Barletta et al. 2019). Although considered less harmful than other anthropic activities (Páez-Osuna 2001), shrimp farming represents a high risk to environmental conservation. The implementation and adverse impacts due to these operations include mangroves devastation, modification of natural drainage, intensifying saltwater intrusion and destruction of aquatic species habitat (Ferreira & Lacerda 2016; Ramos e Silva et al. 2017).

However, the main impact of shrimp farms is caused by discharge of effluents during daily water exchange and harvest – when shrimps are removed after reaching commercial size. These effluents contain high concentrations of nutrients, suspended solids and organic matter (Trott & Alongi 2000; Bui et al. 2013; Herbeck et al. 2013; Cardoso-Mohedano et al. 2016a), causing depletion of dissolved oxygen and increasing water turbidity (Jackson et al. 2003; Barraza-Guardado et al. 2013). Water exchange operations can also modify the morphological balance of coastal systems, especially when it occurs synchronously with tidal oscillation (Roversi et al. 2020). Factors such as location, management, and stocking density (Cardoso-Mohedano et al. 2016b) influence shrimp farm effluent quality and degradation levels. Adjustments in discharges according to estuary residence time and tidal regime can mitigate the impacts of the activity (Páez-Osuna, 2001; Cardoso-Mohedano et al. 2016a).

In general, research about shrimp farming impact on coastal ecosystems focuses on the effluent's composition or the effect of a single farm discharge. Using computational modeling, Cardoso-Mohedano et al. (2018) observed that shrimp farms' discharges induced salinization increasing at least 3.0 psu at a subtropical coastal lagoon for up to 2 tidal cycles. Ferreira et al. (2011) monitored the effluents and pond water on a shrimp farm, observing seasonal variations in water quality, and demonstrated that farm management should be adapted to each season. Ramos e Silva et al. (2010) assessed shrimp farm effluents under intensive, semi-intensive and organic management, observing that variations in organic residue management lead to different water quality standards, more intensive systems lead to worse water quality. Little is known about synergistic effect of multiple ponds discharge. Defining the harvest time series for multiple farms is a challenge, due to the difficulty in acquiring data and quantifying all discharges.

Nutrients from shrimp farming and urban areas are the most common pollution sources in estuaries and coastal lagoons in Northeastern Brazil (Barcellos et al. 2019). The Potengi/Jundiaí Estuary (PJE), on the Northeastern Brazilian coast, is under anthropic discharges that include aquaculture, industrial and domestic wastewater, along surface runoff from agricultural areas. This estuary, where commercial shrimp farming first took place in Brazil, has several shrimp farming ponds along its banks that have promoted marked changes in the ecosystem, which contributed to devastation of at least 30% of the mangrove (Souza & Ramos e Silva 2011). In the present work, the main goal was to analyze the influence of multiple shrimp farms discharges on water quality of the PJE with different strategies of harvest, including a scenario with harvest forecasted by a stochastic model (Roversi 2018). A quantitative index, TRIX, was used to evaluate the level of eutrophication in the PJE for different strategies of harvest, quantifying the trophic conditions of the estuary and indicating the best harvest strategy. Thus, such a water quality model can be employed to support the choice of strategic decisions concerning shrimp farms effluent in coastal areas.

The PJE is located on the coast of Northeastern Brazil (5°43′S to 5°53′S and 35°09′W to 35°21′ W), in Rio Grande do Norte state, integrating the cities of Natal, Macaíba and São Gonçalo do Amarante (Figure 1). The study area presents a tropical climate with dry summers (Alvares et al. 2013), a short annual temperature range and an average temperature of 26.5 °C (INMET 2018). The average annual rainfall is 1,700 mm, with the rainy season between February and June, while average annual evaporation is 2,000 mm (INMET 2018). The upper estuary, near the Potengi river mouth, suffers the greatest environmental impact provoked by anthropic activity (Souza et al. 2016).

Salinity ranges between a minimum of 29.5 psu during the rainy season, and 34.5 psu during the dry season (Miranda et al. 2005). According to vertical salinity stratification patterns, PJE is considered moderately or partially mixed (Frazão 2003). Water temperature pattern is mixed, with an average of 28 °C, typical of low latitude regions (Frazão & Vital 2006). The PJE tributaries are Jundiaí, Potengi and Doce rivers, in addition to other small contributors. These intermittent rivers drain an area of around 5,000 km² (Boski et al. 2015).

Shrimp farms are located along the extension of PJE. Ponds volumes were estimated through the ponds' geometry, considering 1 m depth (Cunha 2010) and surface areas obtained by a RapidEye satellite image taken on 5/31/2014 (MMA 2016). Ponds' areas were compared with surveys conducted in the region by the Institute for Sustainable Development and the Environment (in Portuguese, IDEMA) in 2007 and Google Earth images of 2009, identifying ponds active in 2009. Shrimp farms' harvests were estimated in three 4-hour cycles and the cultivation characteristics were obtained from questionnaires applied in 2009 (Cunha 2010) (Table 1). The ponds are usually emptied by gravity flow through channels that discharge untreated effluents directly into the estuary, starting 2 hours before low tide.

The hydrodynamic and water quality models used in this study are part of the Hydrodynamic Environmental System known as SisBAHIA^® (www.sisbahia.coppe.ufrj.br). The SisBaHiA^® (from Portuguese, Sistema Base de Hidrodinâmica Ambiental) is a system developed by the Coastal and Oceanographic Engineering Department, Oceanic Engineering Program, Federal University of Rio de Janeiro (COPPE/UFRJ), which has been widely used to simulate coastal and estuary ecosystems (e.g. González-Gorbeña et al. 2015; Aguilera et al. 2020; Silva et al. 2020). SisBaHiA^® adopts the finite difference method for time discretization and the finite element method for space discretization. The Water Quality Model (WQM) uses an Eulerian approach, the same spatial grid applied for the hydrodynamics model and allows for the use of different time step lengths in the analyses. The advantage of coupling the two models, as in the current work, appears in the determination of the flow velocities and turbulence coefficients, which was done previously in the hydrodynamic model, that can be used directly in the water quality model.

In WQM, BOD indicates the presence of organic matter and Chlorophyll a is the substitute parameter for phytoplankton biomass. Chlorophyll a is not an expression of production dynamics, but only a general indicator of phytoplankton biomass.

The DO mass balance is influenced by almost all the processes in the system. The oxidation of organic matter consumes DO whereas reaeration process and photosynthetic activity produce it. Other sinks are the nitrification and sediments' oxygen demand. The source of BOD is the death of phytoplankton and HZ. The sinks are oxidation and sedimentation.

Nitrogen and phosphorus are returned to the organic compartment (DON and DOP) via phytoplankton and HZ respiration and death. After mineralization, the organic form is converted in the dissolved inorganic form available for phytoplankton growth. The sources of DON are the death of phytoplankton and HZ. The sinks are mineralization and sedimentation. The sources of NH₃ are mineralization, death of phytoplankton and HZ. The corresponding sinks are nitrification and the growth of phytoplankton. The source of NO₃ is nitrification and the sinks are denitrification and the growth of phytoplankton.

The source of DOP is the death of phytoplankton and HZ and the sinks are, again, mineralization and sedimentation. The sources of DIP are mineralization and death of phytoplankton and HZ and the sinks are precipitation and the growth of phytoplankton.

The source of Chl-a is the growth of phytoplankton, which depends on temperature, availability of nutrients and light, and the sinks are sedimentation, excretion, death and grazing by HZ (herbivory). The herbivorous zooplankton assimilates the ingested phytoplankton (source) and the fraction not assimilated is transferred to the organic matter compartments. The sinks of HZ are precipitation and death.

The modelling area is shown in Figure 1, where the quadratic finite element mesh, which employs foursquare subparametric Lagrangian elements, is included. The bathymetry of the PJE (Figure 1) was obtained from nautical maps from the Directorate for Hydrography and Navigation (DHN 2008, 2009) and from a survey conducted in 2004 by the Geology Department of the Federal University of Rio Grande do Norte (UFRN).

The bottom friction coefficient can be written in terms of the Chèzy coefficient, which depends on the amplitude of equivalent bottom roughness. The amplitude of equivalent bottom roughness was obtained from Frazão & Vital (2006). Sand is predominant in the main channel (approximately 0.030 m), whereas sandbanks predominate near the island and near the Natal Harbor, with a predominance of clayey-silty sediments of around 0.015 m near the river mouths.

The hydrodynamic model was calibrated with data related to the wind and tidal levels from 2009. As the river discharges have little influence on the hydrodynamic flow into PJE, monthly average discharges were used (CAERN 2006) for Potengi and Jundiaí rivers (Table 2), and a permanent flow of 2 m³/s was estimated for Doce river (Ramos e Silva et al. 2006). Domestic and industrial wastewater discharges (Nicodemo 2010) were considered permanent during the simulation (Table 3).

Wind conditions were assumed unsteady but spatially homogeneous. Data used in the model included the time series measured at the Automatic Surface Meteorological Station (INMET 2018), near the estuary, with an hourly measurement. Wind measurements presented the most frequent winds coming from the E-SE direction, with mean velocity equal to 4 m/s, and maximum values of 8.0 m/s. The wind rose representing the time series of the wind from the year 2009 is illustrated in Figure 1.

The PJE hydrodynamic circulation is tidally dominated, with a semi-diurnal tidal regime and maximum tidal amplitude of 2.80 m, characterized as mesotidal (Kumar et al. 2018). In the sea level boundary condition, the free surface position values are prescribed by applying the astronomical tide curve, obtained from 24 harmonic constituents established for the Natal Harbor (FEMAR 2000). The time interval used in the hydrodynamic circulation simulation was equal to 20 seconds, which corresponds to an average Courant number of 3.0. The model for PJE had a warming up period of two tide cycles, to minimize any influence from the initial conditions.

The PJE receives loads of nutrients and organic matter from at least three anthropic origins: urban, industrial and shrimp farms. The WQM was used to analyze the impact of aquaculture effluents in different harvest scenarios. The concentrations of BOD, DO, DOP, DIP, DON, NH₃, NO₃, Chl-a from rivers, sewage treatment plants (STP) and shrimp farms were used as model boundary conditions (Table 3). The concentrations of rivers and STP were obtained from field data collected in 2009 (IDEMA, 2009a, 2009b, 2010; Nicodemo 2010) and chlorophyll a values measured in waste stabilization ponds (Mara et al. 1992). Water quality parameters of shrimp farms, according to the cultivation systems, were entered from field data (Cunha 2010) and literature recommended values (Ferreira et al. 2011; Ma et al. 2013).

Four stations (Figure 1) were defined to analyze WQM results, representing each zone of the estuary: in the upper estuary, S01 – near the Jundiaí river mouth, with a predominance of shrimp farms and ST2 – an intermediate area between the Jundiaí and Potengi rivers mouths, also with a predominance of shrimp farm discharges; in the middle estuary, S03 – near the Potengi river mouth, a region with a strong occurrence of domestic and shrimp farm effluents; and in the lower estuary, S04 – near the estuary mouth, without direct discharges and under the strong influence of the Atlantic Ocean.

The water quality model was simulated for the period of 365 days, with flow velocities and turbulence coefficients defined for the year 2009, previously simulated by the hydrodynamic model. The time step of the WQM was 150 s. Although field data are not sufficient for WQM calibration, water temperature, salinity, DO, BOD, total nitrogen (TN = NO₃ + NH₃ + DON) and total phosphorus (TP = DOP + DIP) concentrations were compared with a few field data collected (IDEMA 2009b, 2009a, 2010; Programa Água Azul 2009, 2010) at several points along the estuary to adjust WQM. The concentrations and flows shown in Table 3 were used to calibrate the water quality model, these values being the most representative of the loads present in the PJE. The adjustments of the constants in the WQM were carried out manually to obtain numerical results close to the measured data. The parameters used in the calibration are those defined by literature (Cunha et al. 2018; Rosman 2018), as there are no specific values for the PJE region.

Four scenarios with varying harvest patterns were developed to evaluate the shrimp farms synergistic impacts and to identify how these effluents changed the water quality of PJE. A reference scenario (RF) was simulated without shrimp farm discharges, with only river loads and sewage treatment plants. Three scenarios adding shrimp farm discharges were simulated: (HG) with harvest flows forecasted by the Q-harvest generator; (ST) with harvests synchronous to the spring tides; (NT) with harvests synchronously to the neap tides. In the scenarios with harvests concentrated during spring (ST) and neap tides (NT), all farms monthly discharge 25% of the ponds. Shrimp farm effluents are also discharged during daily exchanges to maintain oxygen levels in the pond water (Thomas et al. 2010; Aschenbroich et al. 2015).

The mean harvest flow in HG scenario (Figure 2) was 4.13 dam³/h and maximum peaks of 13.5 dam³/h. Harvest flows in the ST and NT scenarios were 116.50 dam³/h. The simulation period made analyses possible under different patterns of meteorological, oceanographic and water quality conditions, including shrimp pond discharges in the entire modeling domain. Differences observed among the different scenarios (RF, HG, ST and NT) on each station were compared by using a two-way analysis of variance (ANOVA). Statistical differentiation between HG-RF and ST-NT scenarios were analyzed by using the Student's t test. The statistical significance was defined as p < 0.05.

The results for free surface position were calibrated by comparison with tidal harmonics prediction data obtained by SisBaHiA^® tide forecast module. Figure 3 present a comparison between the free surface elevation data and that obtained by the numerical simulation. Computed results present agreement with data in terms of phase. The computed tidal amplitudes practically match the tide forecast module amplitudes, but they are slightly shifted up between the 10th to the 15th day. This behavior corresponds to the usual dynamic variation in mean water levels in the inside areas of an estuary. In fact, in such cases, dynamic mean water levels tend to be higher on spring tidal days than on neap tidal days. The accuracy of fit was quantified by the correlation coefficient R² of the regression line between the model results and the field data was computed. The best calibration required an R² value as close as possible to 1.0. The correlation coefficient value (R²) being equal to 0.9893, indicating agreement between simulated and field data. This regression shows a tendency of the model to be less accurate for the minimum level values in the tidal cycle.

The circulation in PJE is characterized by tidal currents, with some influence from fluvial discharges. The spatial pattern of currents shows that the main channel tends to direct the current field in a direction parallel to the contours, increasing their values. At ebb tide, the currents are stronger than at flood tide, due to the influence of lateral and bottom attrition, mainly in the strangled regions. The effect of shallow water is considerable in current variations, being responsible for asymmetries in the distribution of ebb-flood current. The isoline map of occurrence of current velocities below 0.10 m/s for 365 days was extracted to identify possible stagnation areas (Figure 4). Maximum occurrences were found in the upper estuary (70–90%) and in the side arms (80–100%), where water exchange by tidal effects is relatively smaller. The upper estuary is more vulnerable to receiving aquaculture discharges since the low velocities favor nutrient accumulation and can result in long-term eutrophication of water bodies. Minimum occurrence (∼0%) were found in the main channel, in the middle and lower estuary, demonstrating that accumulation of effluent discharges in these areas tends to be lower. However, stagnation areas can possibly develop near the banks, where aquaculture and sewage treatment plant effluents are discharged.

Field data available are not sufficient for a complete calibration, and adjustments were carried out to obtain numerical results close to measured data, mainly for re-aeration and de-oxygenation constants. Figure 5 shows comparisons at station S01 for temperature (T), salinity (S), DO, BOD, TN and TP. The model reproduces temperature variations. Differences between measured and numerically computed salinity can be explained by rivers' daily fluctuations. For DO, BOD, TN and TP, the model reproduced the seasonal variations (rainy and dry). Divergences are certainly originated in boundary conditions, particularly in the determination of loads from domestic, industrial and shrimp farms effluents and from the imposed values of DO concentration from the rivers, assumed constant and no seasonal variations. For a better calibration, a larger temporary series would be necessary.

Table 4 shows a comparison between the results for the different shrimp farm discharge scenarios (HG, ST and NT) and the reference scenario (RF) - without discharges from the shrimp farms. Results indicate that shrimp farms effluents enriched with nutrients and organic matter are accumulated in the upper estuary and side arms. BOD, TN, TP and Chl-a levels were higher in the station near the Jundiaí river mouth, decreased in the station S02 and increased from the middle estuary in all scenarios. The station S04 showed the lowest differences between scenarios with shrimp farm (HG, ST and NT) when compared to RF and was the less impacted by the shrimp farms releases.

All stations showed significative differences (p < 0.05) between HG and RF scenarios. The HG scenario showed higher BOD values in all stations in comparison to the RF scenario. In the upper estuary (S01 and S02), TN, TP and Chl-a concentrations are higher in the HG scenario compared to the RF scenario, corroborating the hypothesis that shrimp farm effluents enriched with nutrients and organic matter are accumulated in the upper estuary. The middle estuary region (S03) presents the highest values for BOD, TN, TP and Chl-a concentrations, indicating that, in this region, the occurrence of low velocities generates stagnation areas, favorable to accumulation of pollutants. Changes in the concentrations of TN and TP may be associated with local water circulation. The presence of effluents from shrimp farms in the HG, ST and NT scenarios imposes an increase in the circulation close to the margins, improving exchanges and reducing stagnation, resulting in slightly water quality.

Lowest BOD was found in the upper estuary without shrimp farm release (RF scenario) at station S02 (2.57 ± 0.85 mg/L). Highest BOD was found in the middle estuary, S03 (6.47 ± 0.69 mg/L), in the HG scenario. The HG and RF scenarios showed few differences in the minimum and maximum BOD values in the S03 station, despite the BOD mean value being 0.52 mg/L higher in the HG scenario.

The most statistically significant differences between the no discharge scenario (RF) and shrimp farms discharge scenarios (HG, ST and NT) occurred in the upper estuary, where velocities are the lowest. The one-way ANOVA indicates significant differences (p < 0.05) between the HG, ST, and NT scenarios in the S01 and S02 stations for all water quality parameters. Higher BOD mean values were found at HG in S01 (5.17 mg/L) and S02 (3.57 mg/L) compared to ST (5.02 and 3.44 mg/L) and NT (5.02 and 3.44 mg/L) scenarios. Total nitrogen also showed more significant differences in the HG scenario in stations S01 (1.74 mg/L) and S02 (1.32 mg/L) compared to ST-NT scenarios (1.71 and 1.29 mg/L). The one-way ANOVA showed no significant differences (p > 0.05) between nutrients and Chl-a in HG scenarios compared to ST and NT scenarios in the middle and lower estuary (S03 and S04). These stations show only BOD significant differences. The HG scenario presented a mean BOD (6.47 mg/L) higher than the ST (6.39 mg/L) and NT (6.40 mg/L) scenarios at S03.

The t-test showed some significant differences between ST and NT scenarios (P < 0.05) only on DBO values at the upper estuary (S01 and S02) and on TP values at station S02. In the S01 station, the ST scenario has a minimum BOD value (3.63 mg/L), which is higher than the value found in the NT scenario (3.58 mg/L), but with equal mean values (5.02 mg/L). The mean BOD value in the S02 station is higher in the ST scenario (3.44 mg/L) compared to the NT scenario (3.42 mg/L). Although being considered statistically significant, these differences do not seem be relevant considering the magnitude of the water quality parameters.

At all stations, the low correlation between TP and Chl-a (R² = 0.69) suggests that the former is not a limiting factor to the growth of phytoplankton. In this case, the limiting factor may be nitrogen, which showed a greater correlation with Chl-a (R² = 0.81). Changes in the concentration of Chl-a may be associated not only with total nitrogen and phosphorus levels but also with water circulation.

BOD concentrations after 30 days of shrimp pond discharges into the estuary during the dry season at station S01 are illustrated in Figure 6. The results show an increase in concentrations immediately after harvests concentrated during the spring (ST) and neap tide (NT) scenarios. The concentrations were reestablished in the water body around 7 days after spring and neap tide discharges. Concentrating discharges during the spring and neap tides should lead to higher concentrations, but this does not in fact occur.

Figure 7 shows the trophic state index (TRIX) results for the HG scenario: average of 5.48 ± 0.90, for all stations; minimum value of 3.01 on 17 February, for station S04; and maximum value of 7.09 on 22 January, for station S03. Station S03 has the highest values, where the eutrophication level was considered high and the water quality was poor, with values greater than 5 for most of the time. In station S04, the PJE waters shows a medium level of eutrophication, with good water quality.

According to the relationship between TRIX scale and ratings of eutrophication level and water quality, the PJE waters shows a medium level of eutrophication, with the water quality ranging from good (TRIX: 4–5) to bad (TRIX: 5–6). TRIX results is consistent with the stations pattern in dissolved nutrients described above, which indicate an accumulation of nutrients in the upper estuary. Daily values of TRIX index at station S02 (Figure 7(b)) show that there is no significant variation (p > 0.05) between the HG (5.18 ± 0.25), NT (5.16 ± 0.26) and ST (5.15 ± 0.26) scenarios. The RF scenario (4.88 ± 0.31) presented TRIX values below 5.0 for 221 days, which represents 61% of modeling time, corresponding to good water quality and medium eutrophication. The values for the TRIX index showed a tendency to rise when penetrating the estuary, especially close to rivers. The lowest values were observed at high tide and the lowest at low tide in all seasons.

Synchronized discharges during spring or neap tide shows a few benefits or differences, but not justifying the restriction of the discharges to any of these periods. These results are inferior to those obtained by Cardoso-Mohedano et al. (2016a), who reported that maximum nutrient concentrations in discharges during the spring tide are up to 10% higher than in the neap tide. The results suggest that nutrient discharges are rapidly assimilated in the main channel and that shrimp farms' effluents increase the chlorophyll a and organic matter concentrations, indicating accumulation of organic matter caused by the activity.

The annual nutrient and organic matter loads for each pollution source were calculated (Table 5) to determine the input of nutrients and organic matter at the estuary. Shrimp ponds released equal organic matter (578.45 ton), total phosphorus (9.63 ton) and total nitrogen (127.23 ton) loads in both ST and NT scenarios, lower than the HG scenario (719.63, 11.96, 158.24 tons, respectively). Comparing to other anthropic pollution sources, sewage treatment plants discharges the highest annual loads into PJE. Lagoa Aerada STP and Baldo STP released more pollutants into the estuary (∼70% of the organic matter and ∼85% of the nutrients) than all shrimp farms together (∼17% of the organic matter and ∼10% of the nutrients). The combination of loads released by shrimp farms is the third highest source of organic matter, TP and TN.

Several studies (Trott & Alongi 2000; Barraza-Guardado et al. 2013) indicate that shrimp farm effluent discharges may not increase nutrient concentrations in the receiving water bodies, but do raise organic matter and chlorophyll a concentrations, especially in the short term and near the farms. By contrast, other studies (Paula Filho et al. 2015; Barcellos et al. 2019) show higher nutrient concentrations than their domestic and industrial effluents. This discrepancy may be caused by differences in local hydrodynamic circulation, the natural variability of nutrients (Cardoso-Mohedano et al. 2016a) and the amount of anthropic effluents.

The PJE showed medium eutrophication and nitrogen was the limiting nutrient. The upper estuary is the most vulnerable area to eutrophication. Shrimp farm discharges increased the level of nutrients, chlorophyll a and dissolved organic matter in the upper estuary. Low current velocities, and small width and depth, please nutrient trapping. Harvests in neap and spring tides resulted in a small improvement in water quality. In these scenarios, the water quality impact remained up to 7 days after the last harvest. The distributed harvests resulted in more harvests during the year, with more annual loads of nutrients and organic matter. Comparing to domestic effluents, loads released by shrimp farms are lower, as the effluents have lower concentrations and flows. The methodology applied may support coastal management, predicting and mitigating the effects of pollution caused by nutrients and organic matter in water bodies where the activity is more relevant. The harvest time series methodology is easy to implement and makes it possible to simulate different effluent discharge scenarios at the end of the production cycle.

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) - Finance Code 001. We thank the National Program for Academic Coordination (PROCAD, project no. 88881.068511/2014–01, Protocol 071/2013), and the project entitled ‘Water quality and circulation in coastal ecosystems: monitoring and modeling for decision-making’.

All relevant data are included in the paper or its Supplementary Information.