Keywords

1 Introduction

The coastal zone of Cartagena, Colombia is inhabited by small coastal communities which have traditionally depended on artisanal fishing for their livelihood. However, these communities are vulnerable to environmental changes that have occurred in Cartagena Bay in terms of water and sediment pollution [1]. These fishing communities have witnessed drastic declines in their fishing resources over the past decades which have led to reduced catch and fish-size [2], while studies have also demonstrated impacts of metal contamination in the fish [3, 4]. As a result, these populations have gradually reduced their diet of fish, though this can have repercussions in terms of health indicators, and they have also shifted their economy towards beach tourism [4]. While tourism may be economically beneficial to the communities now, the sustainability of beach tourism in coastal Cartagena is fragile due to a rapidly increasing number of beach users and limited management infrastructure, while the decline of fishing activities threatens the loss of the traditional culture of artisanal fishing.

There are various issues of water quality in Cartagena Bay that could be affecting the fisheries, including turbid freshwater plumes, eutrophication, hypoxia, bacteria and metal pollution [1], which come from multiple sources of continental runoff, domestic and industrial wastewater [5], making mitigation a complex and challenging task. Such impacts affect fisheries due to the effect of degraded water quality on coral reefs [6], which provide fish with important habitats for feeding, refuge and reproduction. Fish are also impacted directly by insufficient oxygen levels in the water column which limit their respiration [7, 8]. Recent research of oxygen in the water column of Cartagena Bay [1] shows hypoxic conditions as levels of dissolved oxygen and oxygen saturation were found below recommended threshold levels of 4 mg/l [9] and 80% [10], respectively, at depths of just 5–10 m from the surface (Fig. 1). A large input of freshwater in the bay leads to strong vertical stratification which inhibits vertical mixing and thus limits the oxygenation of the water column [11]. Though this problem is alleviated during the dry/windy season (Jan-April), hypoxic conditions are found for most of the year and worst during the transitional season (June-August) when heightened temperatures reduce oxygen levels even further (Fig. 1).

Fig. 1.
figure 1

(see Tosic et al. [1] for details).

Measurements of oxygen profiles in Cartagena Bay during the dry/windy (Feb.), transitional (June) and rainy (Nov.) seasons. Dissolved oxygen concentration (mg/l) is shown as a colour gradient. Oxygen saturation (%) is shown as dashed contour lines. Points (•) represent measurement locations at stations (B1, B3–B8) shown in Fig. 3

Full size image

Monthly monitoring in the bay has found average concentrations of biological oxygen demand (BOD) of 1.15 ± 0.90 mg/l [1], exceeding the maximum threshold value of 1 mg/l established for fishing resources in Cuba [12]. While these conditions are due in part to large amounts of organic matter flowing into the bay from runoff and wastewater discharges [5], a substantial input of nutrients also leads to algal blooms [1, 13, 14]. However, the persistence of turbid plumes during most of the year causes primary production in the bay to be light-limited, rather than nutrient-limited, resulting in a lack of productivity below the pycnocline until the dry season when the plumes subside and algal blooms occur (Fig. 2) [1, 13, 15,16,17]. This seasonal dynamic of primary productivity is also reflected by higher levels of BOD and turbidity found in the bay’s bottom waters during the dry/windy season and the transitional season that follows [1]. These seasonal blooms, in combination with the bay’s stratification (weak vertical mixing) and morphology (small seaward straits, deep bay), result in oxygen depletion in the bottom waters [13, 16, 17].

Fig. 2.
figure 2

(see Tosic et al. [1] for details).

Interpolated maps of surface water turbidity (NTU; above) and bottom water chlorophyll-a (μg/l; below) measured in Cartagena Bay and averaged over the 2014 rainy (left) and 2015 dry/windy (right) seasons

Full size image

To better understand these complex processes leading to hypoxic conditions and evaluate strategies to mitigation the loss of fisheries, numerical models can be used to simulate the hydrodynamic and water quality conditions of the bay. Various studies have been done in Cartagena Bay to model its hydrodynamics and sediment transport [15, 18, 19], eutrophication and oxygen regimes [13, 16,17,18]. Recent modelling work has since applied the MOHID-Water model in Cartagena Bay [11].

Light attenuation in a system such as Cartagena Bay is particularly important for the successful application of a hydrodynamic-water quality model. Given the bay’s characteristics of high surface turbidity and strong vertical stratification, light attenuation plays an important role in the processes of short-wave radiation penetration and the resulting heat flux, which are significant factors for the effective functioning of a hydrodynamic model. Furthermore, the penetration of light is also essential to primary production in such a light-limited system, and thus understanding light attenuation is an important step towards applying a eutrophication model.

In this study, we investigate the processes of light attenuation through monthly field measurements of temperature, salinity, total suspended solids (TSS), turbidity, chlorophyll-a and Secchi depth in Cartagena Bay and apply these data to numerical simulations with the MOHID-Water model. Multiple methods of calculating the short-wave light extinction coefficient (Kd) are carried out utilizing the field measurements and tested in model simulations. By comparing modelled thermoclines to field observations, the best fitting Kd calculation method is identified for further modelling. This first step towards eutrophication modelling in the bay provides an improved understanding of the water quality processes impacting coastal fisheries in Cartagena Bay.

2 Methods

2.1 Study Area

The tropical semi-closed estuary of Cartagena Bay is situated in the southern Caribbean Sea on the north coast of Colombia (10°20’ N, 75°32’ W, Fig. 3). The bay has an average depth of 16 m, a maximum depth of 32 m and a surface area of 84 km2, including a small internal embayment situated to the north. Water exchange with the Caribbean Sea is governed by wind-driven circulation and tidal movement through its two seaward straits: “Bocachica” to the south and “Bocagrande” to the north [11]. Movement through Bocagrande strait is limited by a defensive colonial seawall 2 m below the surface. Bocachica strait consists of a shallow section with depths of 1–3 m, including the Varadero channel, and the Bocachica navigation channel which is 100 m wide and 24 m deep [13, 19]. The tides in the bay have a mixed, mainly diurnal signal with a micro-tidal range of 20–50 cm [20].

Fig. 3.
figure 3

Principal panel: study area showing sampling stations, model calibration points, control points, weather station (SKCG), bathymetry and model domain. Secondary panels: location of Colombia (upper panel); location of the Magdalena River (middle panel); flow of the Magdalena into the Caribbean Sea and along the Dique Canal into Cartagena Bay (lower panel).

Full size image

Estuarine conditions in the bay are generated by the Dique Canal to the south which discharges approximately 50–250 m3/s of freshwater into the bay, the variability of which is strongly related to seasonal runoff from the Magdalena watershed [5, 13]. The flow of freshwater and sediments into the bay generates a highly stratified upper water column with a pronounced pycnocline in the upper 4 m of depth, above which turbid freshwater is restricted from vertical mixing and fine suspended particles tend to remain in the surface layer [11, 13, 16, 17].

2.2 Data Collection

Water quality was monitored monthly in the field from Sept. 2014 to Nov. 2016 between the hours of 9:00–12:00. Measurements were taken from 11 stations (Fig. 3), including one station in the Dique Canal (C0), eight stations in Cartagena Bay (B1–B8) and two stations at the seaward end of Barú peninsula (ZP1–ZP2). At all stations, Secchi depth was measured and CTD casts were deployed using a YSI Castaway measuring salinity and temperature every 30 cm of depth. Grab samples were taken from surface waters while bottom waters (22 m depth) were sampled with a Niskin bottle. Surface samples were collected in triplicate at station C0. A triplicate sample was also taken at a different single station in the bay each month to estimate sample uncertainty. Samples were analyzed at the nearby Cardique Laboratory for total suspended solids (TSS), turbidity and chlorophyll-a by standard methods [21].

At station C0 in the Dique Canal, discharge was measured with a Sontek mini-ADP (1.5 MHz) along a cross-stream transect three times per sampling date from Sept. 2015 to Nov. 2016. Bathymetric data with 0.1 m vertical resolution were digitized from georeferenced nautical maps (#261, 263, 264) published by the Centre for Oceanographic and Hydrographic Research (CIOH-DIMAR). In the 3 × 2 km area of Bocachica strait, the digitized bathymetry was updated with high-precision (1 cm) bathymetric data collected in the field on 17 Nov. 2016.

Hourly METAR data of wind speed, wind direction, air temperature and relative humidity were obtained from station SKCG at Rafael Núñez International Airport (approximately 10 km north of the bay; Fig. 3). Tidal components were obtained for numerous locations offshore of the bay from the finite element tide model FES2004 [22] using the MOHID Studio software. Hourly measurements of water level at a location within the bay were also obtained from the Centre for Oceanographic and Hydrographic Research (CIOH-DIMAR).

2.3 Model Configuration

The hydrodynamics of Cartagena Bay were simulated using the MOHID Water Modelling System [23, 24]. The MOHID Water model is a 3D free surface model with complete thermodynamics and is based on the finite volume approach, assuming hydrostatic balance and the Boussinesq approximation. It also implements a semi-implicit time-step integration scheme and permits combinations of Cartesian and terrain-following sigma coordinates for its vertical discretization [25]. Vertical turbulence is computed by coupling MOHID with the General Ocean Turbulence Model [26].

Model configuration for Cartagena Bay was based on an equally-spaced Cartesian horizontal grid with a resolution of 75 m and a domain area of 196 km2 (Fig. 3). This included an offshore area extending 2.3 km off the bay, though only the results inside the bay are considered within the limits of the monitoring stations used for calibration. A mixed vertical discretization of 22 layers was chosen to reproduce the mixing processes of the highly stratified bay by incorporating a 7-layer sigma domain for the top five meters of depth and a variably-spaced (depth-incrementing) Cartesian domain below that depth. An optimal time step of 20 s was chosen.

The collected data of temperature, salinity, canal discharge, bathymetry, tides, winds and other meteorological data were used to force the hydrodynamic model. This process included a sensitivity analysis to identify the system’s most effective calibration parameters (horizontal viscosity and bottom roughness), and subsequent calibration and validation of the model. For more details on the hydrodynamic model’s configuration, calibration and results, see Tosic et al. [11].

2.4 Evaluation of Kd Calculation Methods

Multiple calculation methods of the short-wave light extinction coefficient (Kd) were compared utilizing the field measurements of total suspended solids (TSS), chlorophyll-a and Secchi depth. These included previously established relationships between Kd-TSS (Eq. 1) and Kd-Secchi (Eq. 7) in Cartagena Bay [15], Kd-TSS and Kd-Secchi relationships in coastal (Eqs. 2, 8) and offshore (Eqs. 3, 9) waters around the United Kingdom [27], a Kd-TSS relationship developed for the Tagus estuary, Portugal (Eq. 4) [28], a relationship for Kd-chlorophyll-a (Eq. 5) in oceanic waters [29], and a combined Kd-TSS-chlorophyll-a relationship (Eq. 6) based on Portela et al. [28] and Parsons et al. [29]. Kd values were computed for all surface stations and months of the monitoring program in order to compare the variability of values obtained using the different calculation methods.

Equation 1.:

Kd = 1.31 * (TSS)0.542

Equation 2.:

Kd = 0.325 + 0.066 * (TSS)

Equation 3.:

Kd = 0.039 + 0.067 * (TSS)

Equation 4.:

Kd = 1.24 + 0.036 * (TSS)

Equation 5.:

Kd = 0.04 + 0.0088 * (Chla) + 0.054 * (Chla)2/3

Equation 6.:

Kd = (0.04 + 0.0088 * (Chla) + 0.054 * (Chla)2/3) * (0.7 + 0.018 * TSS)

Equation 7.:

Kd = 2.3/(Secchi)

Equation 8.:

Kd = e(0.253−1.029*Ln(Secchi))

Equation 9.:

Kd = e(−0.01−0.861*Ln(Secchi))

Hydrodynamic simulations were run using the different Kd values to compare the resulting thermoclines with measurements. Due to limitations of computing time, model simulations completed at the time of writing were limited the use of Eqs. 2, 4 and 5 during the dry/windy season (27 Jan. – 24 Feb. 2016), while continued research will evaluate the remaining Eqs. (1, 3, 6–9) and seasons (rainy, transitional) in the near future. Windy season simulations using Eqs. 2, 4 and 5 were chosen as the first simulations for multiple reasons: The windy season was chosen due to the presence of both suspended sediments and increased chlorophyll-a concentrations during this season, whereas the other seasons tend to have high suspended sediments but low Chla; Eq. 4 was chosen as it has already been included as code within the MOHID model, while Eq. 2 was chosen due to its similarity with Eq. 4; Eq. 5 was chosen to evaluate the importance of chlorophyll-a on Kd, as Chla in the only variable in Eq. 5.

Start and end times for each the dry/windy season simulation coincided with the dates of monthly field sessions. Initial conditions of salinity and temperature were defined by CTD measurements made on the corresponding start date. The Kd values used for simulations were computed using the mean value of the given parameter (TSS, Secchi, Chla) over all surface stations measured at both the start-month and end-month. Kd values were thus held spatio-temporally constant during the simulation period. To avert numerical instabilities, a spin-up period of one day was applied to gradually impose wind stress and open boundary forcings [30]. Outputs from the simulations were compared with field measurements of temperature profiles on the respective end date at six stations points in the bay (Fig. 3). Model performance was quantified by calculating the root mean squared error (RMSE).

3 Results

3.1 Field Observations

Concentrations of total suspended solids (TSS) at the bay’s surface had a range of 2.1–64.8 mg/l with an average of 19.9 ± 13.4 mg/l. Greater concentrations were observed during the months of April-May, Aug. and Oct.–Nov. (Figure 4). Stations B3, B4 and B5 yielded higher concentrations in general due to their proximity to the Dique Canal’s outlet.

Fig. 4.
figure 4

Monthly surface measurements of total suspended solids (TSS), turbidity, Secchi depth and chlorophyll-a at 8 stations in Cartagena Bay (see Fig. 3 for locations).

Full size image
Fig. 5.
figure 5

Average values of the short-wave light extinction coefficient (Kd) calculated from Eqs. 1–9 using observed data of TSS, Chla and Secchi depth. Error bars represent standard deviations.

Full size image

Average concentrations of chlorophyll-a (Chla) of 3.5 ± 4.0 μg/l, and a median value of 2.0 μg/l, were found in the surface waters of Cartagena Bay. A wide range of concentrations of 0.5–22.0 μg/l was observed with peak values of 10–22 μg/l measured between the months of Feb.–May. The occurrence of these peaks during the dry season demonstrates the process of algal blooms occurring following the rainy season when sediment plumes recede and transparency improves. Greater Chla concentrations were found to the north of the bay at stations B4–B8, likely due to improved transparency further from the Dique Canal.

The range of turbidity values found in the bay’s surface waters was 0.7–65.3 NTU with an average of 7.4 ± 9.8 NTU. Greater turbidity levels were found at stations B3, B4 and B5, similar to the spatial pattern of TSS. Turbidity was particularly high in the months of May, June, Aug., and Nov. Increased turbidity during the months of May, Aug. and Nov. reflects the findings of TSS. However, higher turbidity in the month of June occurs after the increases of TSS and Chla found in May.

Measurements with the Secchi disk yielded depths of transparency ranging from 0.1–6.0 m, with an average of 1.8 ± 1.3 m. Greater depths were observed at the stations B7–B8 to the north of the bay, where the lowest turbidity levels were also found. The smallest transparency was observed during the months of May, June, Aug., and Nov., in close agreement with turbidity results. Meanwhile, the peak transparency was observed by Secchi depths in the month of Jan., when the lowest levels of turbidity, TSS and Chla were all found as well.

3.2 Short-Wave Light Extinction Coefficients

The average short-wave light extinction coefficients calculated from the various equations ranged from almost zero to 6.3. The highest coefficient values resulted from Eq. 1 based on TSS with an average Kd value of 6.3 ± 2.1. Also based on TSS, Eqs. 2, 3 and 4, all produced similar average Kd values between 1.4–1.9 with standard deviations between 0.5–0.9. The inclusion of Chla as a variable in Eqs. 5 and 6 resulted in much lower average Kd values of 0.2 ± 0.1.

Calculations of the short-wave light extinction coefficient with Eqs. 7–9 using measurements of Secchi depth resulted in high variability as standard deviations were greater than the average values themselves. An average coefficient of 2.7 ± 4.5 was calculated with Eq. 7, while lower Kd values of 1.6 ± 2.7 and 1.1 ± 1.3 were found with Eqs. 8 and 9, respectively.

3.3 Modelling Results

Results of the model simulations yielded relatively good predictions of temperatures when compared to observations, with RMSE values of 0.33 ℃, 0.35 ℃ and 0.37 ℃ for simulations using Eqs. 2, 4 and 5, respectively. Results of the model simulation using Eq. 5 for calculation of the Kd coefficient based on concentrations of chlorophyll-a were found to be least effective in reproducing the temperature profiles of the dry/windy season (RMSE = 0.37 ℃; Fig. 6). Temperature profiles generated by this simulation at different stations in Cartagena Bay were less vertically stratified than measurements. Within the top meter of surface water, the predicted temperature profiles were vertically uniform at stations B4–B7 suggesting a lack of light attenuation at the surface, in sharp contrast to the observed temperature profiles. At stations B1 and B3, the temperature profiles predicted by the simulation using Eq. 5 were inverse in the top 4 and 2 m of water, respectively, showing a drastic dissimilarity to the observations. However, model results for the subsurface seemed to improve significantly using Eq. 5, which could reflect an agreement between this Chla-based equation and the observations of increased Chla in subsurface waters during the dry/windy season (Fig. 2).

Fig. 6.
figure 6

Model results (Sim) using different Kd equations (Eqs. 2, 4, 5) compared to measurements (Obs) of temperature profiles at the end of the dry/windy season simulation at six stations in Cartagena Bay (B1, B3–B7).

Full size image

Predicted results of temperature from the simulations using Eqs. 2 (RMSE = 0.33 ℃) and 4 (RMSE = 0.35 ℃) were more similar to observations than those of Eq. 5. With Eqs. 2 and 4, the predicted temperature profiles were more similar to observations in terms of the pattern of vertical stratification, with the exception of predicted temperatures at station B1 where an inverse pattern was also found at depths of 2–3 m in the simulated results. However, in all cases there is a difference of approximately 1 ℃ between in the observed and predicted temperatures near the surface, showing that the model is underestimating surface water temperatures. This discrepancy is reduced gradually to approximately 0.1–0.2 ℃ as depths increase. The simulation using Eq. 2 performed slightly better than that of Eq. 4 in terms of reproducing observed temperature profiles, though the difference in performance is quite small.

4 Discussion

4.1 Comparison of Kd Calculation Methods

The improved performance found using Eqs. 2 and 4 when compared to Eq. 5 may be expected given the high concentrations of TSS found in the bay. As Eqs. 2 and 4 include TSS as a variable, and Eq. 5 does not, it would be reasonable to expect that the model would be more dependent on concentrations of TSS than those of Chla, which is the variable in Eq. 5. Perhaps Chla might be a more important variable in simulations run between the months of Feb.–May when higher Chla concentrations were observed. In this regard, additional simulations should be run during the months of Feb.–May using Eq. 5 based on the Kd-Chla relationship or Eq. 6 based on a combination of Kd-Chla and Kd-TSS relationships.

Additional simulations should also be run using Eq. 1, as it was previously developed for Cartagena Bay. Initial simulations were focused on using Eqs. 2 and 4 because of the similarity in average Kd values that they produced (Fig. 5) and given that Eq. 4 is integrated in the MOHID Studio software. However, considering the discrepancy found between observed and predicted temperatures, perhaps the higher Kd values produced by Eq. 1 will produce results that are more consistent with the steep vertical temperature gradients observed in the bay’s surface layers.

The use of Eqs. 7–9 may be expected to produce results similar to Eq. 2 and 4 due to the similarity in the average Kd values calculated by these equations (Fig. 5). The large standard deviations calculated using Eqs. 7–9 show the high variability in Kd values generated with these equations. Therefore, it may be illogical to use an average Kd value calculated with such equations that result in such high variability. In this regard, a more realistic approach with all of the equations would be to use a model configuration that permits a continuous calculation of Kd in space and time (for each cell and at each time-step) according to the variability of the parameters TSS, Chla or Secchi depth.

Limitations in the model’s prediction of the thermocline could also be due to inaccuracies in predicting the river plume with the model’s present configuration. As the freshwater discharge is based on interpolation of monthly measurements and the winds are taken hourly from a single location, these parameters forcing the model could generate inconsistencies with the field observations. These issues could possibly even be exacerbated using a model configuration that permits a continuous calculation of Kd in space and time, as the model’s accuracy in predicting surface water variables such as TSS or Chla (which are strongly dependent on plume dispersion patterns) would then play an even greater role in determining thermocline predictions.

One variable that could possibly contribute to the underestimation of surface temperature of the model could be additional freshwater input not considered in the model’s present configuration. Additional sources of freshwater not considered include a series of smaller canals flowing from the industrial sector along the bay’s east coast, surface runoff from the small catchment area of the bay itself, and occasional discharges from the city’s sewerage system through an outdated submarine outfall and backup outlets along the coast when the system overflows. As this freshwater is warmer than seawater but not accounted for in the model, it may result in a bias for the model to underestimate surface temperatures in the bay,

Additional simulations should also be run for the rainy season and traditional seasons when conditions in turbidity differ from the dry/windy season. Ideally a Kd-relationship that satisfies the varying conditions of the different seasons could be found. However, it may also be possible that the strong seasonal variability in water quality parameters like turbidity, TSS and Chla could result in a seasonal variability in light attenuation processes that require different Kd calculation methods depending on the season.

4.2 Social Implications

Though a better understanding of light attenuation properties in Cartagena Bay may not result in an immediate improvement in the livelihood of artisanal fishermen, an estuarine water model would provide environmental authorities with a powerful tool for managing the coastal zone. Mitigation strategies are needed to improve the coastal zone’s water quality which is essential for fishing resources. Such strategies should be evaluated with model simulations prior to implementation in order to assess the consequences on coastal water quality, as they may be contrary to the intended result [11]. Given the importance of light attenuation on the processes of primary productivity and oxygen in the bay, this research provides valuable knowledge for eutrophication modeling which would be beneficial to mitigation planning.

However, it should also be noted that in addition to mitigating water pollution, there is also a need for improved fishery management. Research done by Escobar et al. [2] on fish populations in the study area showed that genetic diversity in fish populations was high, which supports the sustainability of the resource. However, the authors also found that the mean body length for all species was significantly smaller than body length at maturity. While this latter finding implies that fish stocks in this coastal zone are low, it also indicates a need for improved fishery management, as fishermen should not be consciously capturing immature fish. Given the recent findings of Tosic et al. [1, 5] and those of Escobar et al. [2], the cause of this low abundance in the fish stock is likely due to both pollution and overfishing.

5 Conclusions

Estuarine modelling in a semi-enclosed water body such as Cartagena Bay requires an improved understanding of processes of light attenuation. The bay is characterized by turbid surface waters and hypoxic conditions throughout most of the year. However, during the dry/windy season (Jan.–April) when greater vertical mixing improves oxygen conditions and receding sediment plumes result in increased transparency and primary productivity in the water column.

Calculation of the short-wave light extinction coefficient can vary greatly depending on the equation used. Average Kd values calculated for Cartagena Bay from nine different equations ranged in values from almost zero to 6.3. Numerical modelling simulations showed that the selection of the appropriate method for Kd calculation is important for the model’s functionality. Among the three equations evaluated, the two equations that included TSS as a variable generated temperature profiles that were more similar to field observations than the profiles simulated using an equation including Chla as a variable. However, further research is needed to evaluate the remaining Eqs. (1, 3, 6–9) and seasons (rainy, transitional).

An estuarine water model presents environmental authorities with a powerful tool for managing the coastal zone as it allows for the evaluation of mitigation strategies needed to improve coastal zone water quality. Such strategies are particularly important to the artisanal fisheries of Cartagena Bay and the vulnerable coastal communities that depend on them. As the economies, traditional culture and public health of these societies may be impacted by water contamination, the use of water quality models by coastal management in Cartagena Bay would be beneficial to the effective planning of pollution mitigation.