Land Use and Land Cover Shape River Water Quality at a Continental Caribbean Land-Ocean Interface

Clark, Kasey E.; Bravo, Viviana D.; Giddings, Sarah N.; Davis, Kristen A.; Pawlak, Geno; Torres, Mark A.; Adelson, Anne E.; César-Ávila, Carolina I.; Boza, Ximena; Collin, Rachel

doi:10.3389/frwa.2022.737920

ORIGINAL RESEARCH article

Front. Water, 11 March 2022

Sec. Environmental Water Quality

Volume 4 - 2022 | https://doi.org/10.3389/frwa.2022.737920

This article is part of the Research Topic The UN International Day of Forests 2022: Forests and Water View all 4 articles

Land Use and Land Cover Shape River Water Quality at a Continental Caribbean Land-Ocean Interface

$\nKasey E. Clark &#x;$ Kasey E. Clark¹^*^†

Geno Pawlak^2,5

Carolina I. César-Ávila¹

Ximena Boza¹

Rachel Collin¹

¹Smithsonian Tropical Research Institute, Balboa Ancón, Panama
²Scripps Institution of Oceanography, University of California, San Diego, La Jolla, CA, United States
³Department of Civil and Environmental Engineering, University of California, Irvine, Irvine, CA, United States
⁴Department of Earth System Science, University of California, Irvine, Irvine, CA, United States
⁵Department of Mechanical and Aerospace Engineering, University of California, San Diego, La Jolla, CA, United States
⁶Department of Earth, Environmental, and Planetary Sciences, Rice University, Houston, TX, United States

Land use and land cover (LULC) can significantly alter river water, which can in turn have important impacts on downstream coastal ecosystems by delivering nutrients that promote marine eutrophication and hypoxia. Well-documented in temperate systems, less is known about the way land cover relates to water quality in low-lying coastal zones in the tropics. Here we evaluate the catchment LULC and the physical and chemical characteristics of six rivers that contribute flow into a seasonally hypoxic tropical bay in Bocas del Toro, Panama. From July 2019 to March 2020, we routinely surveyed eight physical and chemical characteristics (temperature, specific conductivity, salinity, pH, dissolved oxygen (DO), nitrate and nitrite, ammonium, and phosphate). Our goals were to determine how these physical and chemical characteristics of the rivers reflect the LULC, to compare the water quality of the focal rivers to rivers across Panama, and to discuss the potential impacts of river discharge in the Bay. Overall, we found that the six focal rivers have significantly different river water characteristics that can be linked to catchment LULC and that water quality of rivers 10 s of kilometers apart could differ drastically. Two focal catchments dominated by pristine peat swamp vegetation in San San Pond Sak, showed characteristics typical of blackwater rivers, with low pH, dissolved oxygen, and nutrients. The remaining four catchments were largely mountainous with >50% forest cover. In these rivers, variation in nutrient concentrations were associated with percent urbanization. Comparisons across Panamanian rivers covered in a national survey to our focal rivers shows that saltwater intrusions and low DO of coastal swamp rivers may result in their classification by a standardized water quality index as having slightly contaminated water quality, despite this being their natural state. Examination of deforestation over the last 20 years, show that changes were <10% in the focal catchments, were larger in the small mountainous catchments and suggest that in the past 20 years the physical and chemical characteristics of river water that contributes to Almirante Bay may have shifted slightly in response to these moderate land use changes. (See supplementary information for Spanish-language abstract).

Introduction

Global environmental change, such as anthropogenic and climate changes, impact river water quality and biogeochemical cycles (Carpenter et al., 1992; Schlesinger et al., 2006; Stevenson and Sabater, 2010; Muehlbauer et al., 2020). Changes in tropical climates and ecologies may differ from changes experienced in non-tropical regions (Vitousek, 1994; Keys and McConnell, 2005; Intergovernmental Panel on Climate Change, 2014). The humid tropics cover 1/5 of the global land surface and rivers in this region generate the greatest fraction of global runoff (Fekete et al., 2002). Despite this, however, most of our understanding of rivers comes from research carried out in temperate regions (Turner and Rabalais, 1991; Rosgen, 1994; Smith et al., 1997; Brogna et al., 2017).

Oligotrophic tropical waters support some of the most diverse and productive marine ecosystems such as coral reefs and seagrass meadows. Oligotrophic ecosystems thrive under low nutrient and low turbidity conditions, which promote benthic diversity, complex food webs, environmental stability, and ecosystem spatial heterogeneity (Hallock, 1988). These ecosystems, however, are under threat from global environmental change (Waycott et al., 2009; Poloczanska et al., 2013; Chollett et al., 2017; Sippo et al., 2018), and oligotrophic systems are particularly threatened by nutrient-laden runoff. Under normal conditions, rivers are an important source of vital elements to these marine biogeochemical cycles (Eyre, 1994; Poulton and Raiswell, 2002; Morel and Price, 2003; Bruland and Lohan, 2006; Howarth et al., 2011; Fernandes et al., 2018). However, excessive inputs of nutrients from terrestrial erosion, and agricultural and urban runoff can cause ecosystem degradation, including eutrophication and hypoxia in downstream estuaries and coastal zones (Smith, 2006; Rabalais et al., 2009; Garnier et al., 2010; Alves et al., 2013; Silva et al., 2013; Salas et al., 2015; Bailey et al., 2020).

River water quality is related to the catchment land use and land cover (LULC) (Ngoye and Machiwa, 2004; Slaughter and Mantel, 2013; Kebede et al., 2014; Gallo et al., 2015; de Oliveira et al., 2016; Dumago et al., 2018; Łaszewski et al., 2021). The impacts of LULC change on river water quality are enhanced in the tropics due to biotic factors, such as higher biomass, and more productive tropical forests (Saatchi et al., 2011; Malhi, 2012); a abiotic factors, such as higher precipitation, intense and frequent flooding, and warmer temperatures (Stallard and Murphy, 2012); and anthropogenic factors such as rapidly expanding agricultural practices and population growth which can lead to different water quality responses in tropical rivers compared to temperate rivers (Laurance et al., 2014; Tanaka et al., 2021). LULC changes in the tropics are driven by logging, pasture/ranching, urbanization, agriculture, and burning of vegetation (Tanaka et al., 2021). Generally, LULC change in the tropics increases nitrates, phosphates ( ${PO}_{4}^{-}$ ), ammonium ( ${NH}_{4}^{+}$ ), electrical conductivity (including Na⁺, Mg²⁺, Cl⁻, K⁺, Ca²⁺), chlorophyll-a and suspended solids/turbidity and decreases dissolved oxygen (DO) in downstream rivers and estuaries (Deegan et al., 2011; Riskin et al., 2017; Santos and De Paula, 2019; Tanaka et al., 2021). In some situations, opposite trends have been found. For example river nitrate concentrations were found to decrease in Panama and Brazil where primary forest was converted to pasture (Deegan et al., 2011; Valiela et al., 2013, 2014). In contrast, some river water characteristics such as pH do not show a consistent trend with LULC change (Ling et al., 2016; Ríos-Villamizar et al., 2017; Santos and De Paula, 2019; Tanaka et al., 2021), and underlying geology may be a more important factor in these cases (Young et al., 2005).

Here we describe river water physical and chemical characteristics, including nutrient concentrations of six focal rivers that provide freshwater inputs to Almirante Bay on the Caribbean coast of Bocas del Toro, Panama. This semi-enclosed bay experiences well-documented periodic hypoxia associated with warming that has negatively impacted the local marine life (Altieri et al., 2017; Lucey N. et al., 2020; Lucey N. M. et al., 2020; Weinstock and Collin, 2021). Recent analyses have demonstrated that stratification and temperature inversions associated with hypoxia in Almirante Bay (Neal et al., 2014) are likely driven by physical processes (Adelson et al., in review). However, it is not yet known if or how nutrient inputs might contribute to hypoxia in the Bay. As a step toward understanding how land cover and riverine inputs impact Almirante Bay, we addressed the following questions: (1) How temporally and spatially variable are the physical and chemical (temperature, specific conductivity, salinity, pH, DO, nitrate and nitrites as total oxidized nitrogen (TON), ammonium, and phosphate) characteristics of six focal rivers that contribute to Almirante Bay? (2) How does LULC influence the physical and chemical properties of the focal rivers? (3) How do the focal rivers compare to rivers across Panama?

Materials and Methods

Study Sites

Within the Almirante Bay catchment, there are many small catchments that drain into the Bay, with a combined area of 445 km² (Figure 1). The mountainous terrain to the south of the Bay is mostly secondary forest, with some pasture and scrub/bushes (Figure 2). The low-lying coastal areas are swampy, and an expanse of peat swamp forest lies to the northwest. The peat swamp forest is in a national protected area, San San Pond Sak, an internationally recognized Ramsar wetland (Tejada, 1993; Ramsar, 2020). Ramsar wetlands are wetlands of international importance especially as waterfowl habitat (Ramsar, 2016). Within the Almirante Bay catchment there are three towns, Almirante located on the coast of the mainland in the back of the bay, Bocas del Toro located on Colon Island, and Old Bank, located on Bastimentos Island (Figure 1). We monitored six focal rivers that drain into or near Almirante Bay that were accessible by road or boat, and which span a range of catchment areas and vegetation types (Figures 1, 2). Monitoring sites were selected in each of these rivers in July 2019 (Figure 1). Measurements were taken at the focal rivers every 2 weeks over 9 months from July 2019 to March 2020. Sampling ended abruptly at the end of March, 2020 due to the COVID-19 pandemic (see Supplementary Table 1). The river monitoring took place over 2 days; the small mountainous rivers were measured on land 1 day, and the coastal rivers were measured by boat the subsequent day. Some additional sampling took place at the start of the field campaign in July 2019, and in November 2019 to capture storm flow (see Supplementary Table 1). Based on geomorphological and vegetation differences, we divided the focal river catchments into three general river types. They are as follows: two swamp rivers (Black and Changuinola Canal); one large mountainous river (the Changuinola River); and three small mountainous rivers (Nigua, Oeste, and Pastores Rivers).

FIGURE 1

Figure 1. Almirante Bay and the Bocas del Toro region of Panama. The monitored focal catchments that flow directly into Almirante Bay are outlined in black with red fill and the monitored focal catchment Changuinola flows into the ocean is outlined and filled with yellow. The locations of the river monitoring sites are indicated by filled circles. Other elements of importance are urban areas (Ministerio de Ambiente, 2012), primary roads modified from STRI GIS Laboratory (2011), and the Changuinola Dam. The basemap is hillshade derived from SRTM DEM (CGIAR and NASA, 2019). The inset shows the extent of the Changuinola catchment and the location of the enlarged area is delineated by the black box.

FIGURE 2

Figure 2. Land use and land cover (LULC) for the Almirante Bay and Bocas del Toro region, modified from Ministerio de Ambiente (2012) and Sjögersten et al. (2020). The focal catchments are outlined in dark red, and the San San Pond Sak Ramsar wetland and national protected area (Ramsar, 2020) is represented by thatched gray lines (Tejada, 1993). The inset shows the LULC of the Changuinola catchment, the extent of the Almirante Bay catchment, and location of the enlarged area which is delineated by the black box.

Changuinola River

The Changuinola River is the third largest river in Panama (subsequently referred to as Changuinola) (Comité de Alto Nivel de Seguridad Hídrica, 2016). The catchment boundary runs along the Continental Divide, and stretches into Costa Rica to the west (Figure 1 inset). The river discharges into the open coast of the Caribbean Sea, northwest of Almirante Bay (Figure 1). Under certain conditions, river discharge moves southeast along the shore and enters the bay. The catchment has an area of 3,124 km² and a mean elevation of 1,283 m above sea level (masl) (see Supplementary Table 2). The mountainous portions are mostly mature forest (Figure 2). In the lower reaches and in the floodplain, there is pasture, banana plantations, and the town of Changuinola (Figure 2). There is a large hydroelectric dam in the mountains 46 km upstream draining ~45% of the Changuinola catchment, with a reservoir that extends 15 km upstream (Figure 1). The dam was finished in 2011, it appears to moderate the river flow, generally reducing high and low discharges (ANAM, 2009, 2013). The Changuinola monitoring site was located in the floodplain, ~3 km from the river mouth. The site is influenced by the tide and is located near sea level. The river is wide, shallow, and fast flowing; thus, the river was sampled from an anchored boat (Figure 3a).

FIGURE 3

Figure 3. Focal river monitoring sites at baseflow facing upstream at (a) Changuinola River, (b) Changuinola Canal, (c) Black River, (d) Nigua River, (e) Oeste River, and (f) Pastores River. Photos taken by: V. Bravo (a,c), X. Boza (b,d,f), and C. César-Ávila (e).

Canal

In 1903 the Changuinola Canal (subsequently referred to as the Canal to avoid confusion with the Changuinola) was trenched parallel to the coast on the outer region of San San Pond Sak wetland to move bananas from the Changuinola plantations to Almirante Bay for export (Cramer, 2013). The monitoring site was near sea-level and located roughly halfway from either end of the canal (Figure 1; see Supplementary Table 2). The Canal catchment is dominated by swamp vegetation, with raphia palm and mixed flooded forest (Sjögersten et al., 2020) and is part of the San San Pond Sak wetland (Figure 2). We observed that the flow direction of the canal was variable and likely influenced by the tides, wind direction, and discharge from the Changuinola. Although the bulk of the Changuinola discharges into the ocean, some river water enters the Canal in the Changuinola floodplain (Figure 1). As such, the true contribution area is unknown, but the Canal has a minimum contribution area of 21.8 km² (see Supplementary Table 2). To maintain position while sampling the river water, a rope was tied taut across the river, the motor was turned off, and the boat was kept in place using the rope (Figure 3b).

Black River

The “black” river appears to be unnamed on maps, but for ease of communication we refer to it as the Black River due to its dark color (subsequently referred to as Black). This river is distinct from the nearby Rio Negro, located west of the Changuinola floodplain (ANAM, 2013). The river monitoring site was near sea-level, 2.5 km from the river mouth and has a contributing area of ~46 km² (see Supplementary Table 2). The Black catchment is part of San San Pond Sak wetland and drains into Almirante Bay (Figure 2). It has peat deposits up to 8 m thick and 4,000 years old (Phillips et al., 1997). The majority of this catchment is swampy, with flooded vegetation (stunted forest and sawgrass bog plain), mixed flooded forest, and raphia palm swamp (Troxler et al., 2012; Sjögersten et al., 2020). While most of this catchment is pristine, a small portion to the south west outside of the protected area has been converted to pasture (Figure 2). Despite our efforts to select a sampling site beyond estuarine influence, we found that during the dry season a saltwater wedge reached the site. Under non-storm conditions, the river velocity was very slow and often followed the ebb and flow of the tides (pers. obs.). The river was sampled following the same method used in the Canal (Figure 3c).

Nigua River

The Nigua River monitoring site was situated at ~25 masl, 2.2 km from the mouth, with a catchment contribution area of 13.4 km² and a mean elevation of 128 masl (see Supplementary Table 2). The monitoring site was located under a small bridge accessible by a dirt road (Figure 3d). Most of the catchment land cover is secondary forest, with some scrub/bushes and pasture (Ministerio de Ambiente, 2012) (Figure 2). The monitoring site was upstream of the town of Almirante, with a population of 14,000 (Comité de Alto Nivel de Seguridad Hídrica, 2016), with most of the town located in the low-lying areas (Figure 1). There is a paved road in the catchment that runs alongside the river (Figure 1). There is a drinking water treatment facility with a processing capacity of 7,000 m³ d⁻¹ (Comité de Alto Nivel de Seguridad Hídrica, 2016), ~600 m upstream of the river monitoring site.

Oeste River

The Oeste River monitoring site was situated at ~45 masl, 3.2 km from the mouth, has a contribution area of 14.7 km² and mean elevation of 260 masl (Supplementary Table 2). The land cover in this small mountainous catchment is dominated by secondary forest (Figure 2). The river monitoring site was located under a tall bridge and was accessible by truck. River monitoring was carried out in the river at low flow and from the bridge at high flow (Figure 3e).

Pastores River

The Pastores River monitoring site was situated at ~35 masl and is ~1.6 km from the mouth (Supplementary Table 2). It is the smallest study catchment, with a contribution area of 6.3 km² and a mean elevation of 248 masl (see Supplementary Table 2). The Pastores catchment land cover is dominated by secondary forest and scrub/bushes (Figure 2). There are mangrove forests along the coast (Figure 2). The river monitoring site was located near the intersection of the coastal road (Figure 1). River monitoring was carried out in the river on foot at low to medium flow, and from one riverbank at high flow (Figure 3f).

River Water Characteristics

Physical and Chemical Measurements

We took measurements at 50% of the river depth, in the middle of the river using a YSI Professional Pro Plus Multiparameter Instrument, to measure the following: temperature (°C); specific conductivity (mS/cm); salinity (ppt); pH; and dissolved oxygen (DO) (mg/L). Specific conductivity is often reported in published descriptions of river water (where salinity is not generally reported), and we report it here to aid comparisons with other studies. The YSI was calibrated following the protocols outlined in the manual (YSI Inc., 2009). The conductivity/temperature sensor (YSI 5560 Conductivity, Temperature Probe) was calibrated every 2 weeks for freshwater (YSI 3161 Conductivity Calibrator Solution, 1,000 μS/cm). The pH sensor (Model 1001 pH sensor for ProSeries) was calibrated once every 2 weeks with a neutral buffer (LabChem, buffer solution pH 7.0) and with a basic buffer (LabChem, buffer solution pH 10.0). Once a month the DO sensor (Model 2002 Galvanic DO Sensor, Pro Series) was cleaned to remove moisture and to replace the membrane (YSI 5913 membrane kit 1.25 mil polyethylene) and electrolyte solution (YSI Galvanic Oxygen Sensor, NaCl Electrolyte). The instrument was turned on at least 10 min prior to field measurements to ensure that the DO sensor was acclimatized. The YSI sensor was fixed to a 3 m pole or to a weighted rope and lowered into the water to take the measurement at the appropriate depth.

Nutrients

We monitored phosphates, ammonium, and nitrate + nitrite (total oxidized nitrogen; TON) concentrations. River water was collected at deeper sites with a van Dorn bottle at 50% depth along the vertical profile from the middle of the river (see Supplementary Table 1). A manual grab sample using a 500 ml collection bottle (Nalgene, HDPE) was collected at the shallower sites (<25 cm deep) (see Supplementary Table 1). In the preliminary month of the study, water samples were collected in the field, filtered in the laboratory with GF/F glass fiber membranes (25 mm, 0.7 μm, Whatman, Sigma-Aldrich; see Supplementary Table 1), stored in 60 ml Nalgene bottles, and kept frozen until analysis. For the remainder of the study, river water was syringe filtered in the field directly after collection, with a 0.45 or 0.22 μm (25 mm, Sterlitech) polyethersulfone (PES) membrane filter (see Supplementary Table 1). A new sterile 60 mL syringe (Henke Sass Wolk, silicone oil free) was rinsed three times with the sample. The syringe was filled with sample, the filter was attached and expressed with ~2 ml of the sample to rinse the syringe. Multiple filters were used for samples with high suspended load, where each new filter was rinsed before use. A 50 mL sample vial (Corning, polypropylene centrifuge tube) was rinsed out three times using approximately 3 ml of filtrate per rinse. The sample vial was filled two thirds with filtrate and stored on ice, transferred to a freezer, and stored frozen until analysis.

TON and ammonium were determined through flow injection analysis (FIA) on a Lachat Quikchem 8500 system (Hach Ltd, Loveland, CO) (White, 2003; Tucker, 2007, 2008) in the soils laboratory at Smithsonian Tropical Research Institute (STRI), Panama City, Panama. Glassware and Lachat tubes were cleaned, and acid washed with 10% HCl, rinsed twice with a small amount of defrosted sample which was discarded before adding 10 ml of sample to the tube. The detection limit for TON ranged from 0.005 to 0.01 mg N/L. When TON or ammonium river samples fell below the detection limit, we assigned the sample half of the value of the detection limit for statistical analysis. Filtered laboratory water samples from the Bocas station were run alongside the river samples as controls. The filtered laboratory water samples all had results below the detection threshold.

Orthophosphates ( ${PO}_{4}^{3 -}$ ) were determined through FIA, as described above until October 2019 (see Supplementary Table 1). The detection limit was 0.002 mg P/L, and none of the river samples fell below this limit. Because this approach can be subject to salt error or silicon interference (Murphy and Riley, 1962; Patey et al., 2008; Stetson et al., 2019), we changed from FIA to spectrophotometry when we observed occasional evidence of salt water interference in the brackish river water samples. Thus, from mid-October 2019 onwards phosphates were determined using a molybdenum blue method with absorbance measured on a spectrometer (Hach, DR 5000 Spectrophotometer) (Murphy and Riley, 1962) (see Supplementary Table 1) with a 5 cm quartz cell that was rinsed with sample before use. None of the river samples fell below the detection limit for phosphates, 0.0005 mg P/L. Some of the filtered laboratory water samples run alongside the river samples had values above the detection limit ranging from 0.003 to 0.008 mg P/L, an order of magnitude below the values recorded for the river water samples.

Statistical Analyses

Statistical analyses were carried out on river samples that had been characterized for all eight physical and chemical parameters (see Supplementary Table 1). One-way ANOVAs were applied to each variable independently to determine if the temperature, salinity, conductivity, pH, DO, ${NH}_{4}^{+}$ , TON, and ${PO}_{4}^{-}$ differed significantly among the rivers. When the models were significant, post-hoc Tukey tests were applied to determine which rivers differed significantly from the others. Salinity and conductivity had highly skewed distributions and were therefore log-transformed prior to analysis so they would conform to the assumptions on parametric analyses. All statistical analyses were conducted in JMP 15.0.0.

To determine if temporal changes in water parameters were similar across the different catchments, we calculated the pairwise correlation between each pair of rivers for each variable for each of the fortnightly measurements. There were between 9 and 19 measurements per site that were used in the statistical analysis (see Supplementary Table 1). This approach is not corrected for potential temporal autocorrelation, as we had too few measurements to conduct time series analysis. We applied a Bonferroni correction to insure there was a multiple comparison correction of the p-values. In addition, preliminary data from rainfall gauges and preliminary continuous river level data indicated that, apart from the Changuinola, the flushing time from storm events was <2 weeks. This suggests that our fortnightly sampling rate should result in independent samples. Conductivity measurements were analyzed with non-parametric correlation coefficients due to the un-correctable skewed distribution of values resulting from saltwater intrusion during a few sampling days.

To determine if water from the six focal rivers can be distinguished statistically from each other based on the physical and chemical parameters, we applied a discriminant analysis (DA) using the quadratic difference covariance option in JMP, as some of our data failed to meet the homoscedasticity assumptions (Büyüköztürk and Çokluk-Bökeoglu, 2008). DA is a multivariate approach used to determine what linear combinations of variables most effectively sort observations into predefined groups (i.e., rivers). As we had numerous measurements of only six rivers, this approach focuses on the factors that explain the differences between the rivers (Lachenbruch and Goldstein, 1979). We also conducted a principal component analysis (PCA) to describe the overall patterns of correlation in physical and chemical parameters observed across the entire dataset, and used a PCA conducted on river means to compare with LULC (see below). Analyses were conducted on data from site-week combinations that had no missing data. This had the effect of eliminating the sampling dates where we detected particularly brackish or marine water at the Black site (see Supplementary Table 1).

Catchment and River Delineation

In ArcMap (v 10.8) the hydrology toolset was used to delineate the catchment boundaries and the rivers in the study area (ESRI, 2016) based on the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) (CGIAR and NASA, 2019). The catchment boundaries were determined automatically using the following functions: Fill, Flow Direction, Flow Accumulation, and Basin. Then, the basins were converted from raster to polygons (ESRI, 2016). The Almirante Bay catchment (ABC) contained 1,278 polygons that were manually merged using the editing tool to 14 (Figure 1). The bulk of these automatically generated catchment polygons were very small and/or lined the waterfront. We merged the polygons using the editing feature grouping sub-catchments together based on similar outflow locations into Almirante Bay (Figure 1). Finally, the six focal catchment polygons were exported from the dataset and manually edited based on topography to include only the contribution area of the river monitoring sites (Figure 1).

Rivers were determined automatically from the Flow Direction output following a standard protocol (ESRI, 2016), where rivers were considered to occur where flow accumulation was >100 for 90 m by 90 m cells. The output was converted from raster to polyline (ESRI, 2016) and clipped for the Almirante Bay and Changuinola catchments. Using the editing feature, the river polygons were modified manually, so that Almirante Bay catchment rivers had only one to four rivers per catchment, with many of the small polylines along the coast removed (Figure 1). Some rivers were manually edited or added because the SRTM DEM had low resolution and did not pick up some of the low-lying flat coastal areas, such as Black and Canal. Using the editing feature, these river polylines were edited using our GPS field tracks, and/or using the high-resolution imagery basemap provided in ArcMap (ESRI, 2020) as a visual guide.

Land Use and Land Cover

The land cover map produced for this study consists of two main datasets. The LULC for areas outside of San San Pond Sak were extracted from a country-wide digital LULC map from 2012 produced by the Ministry of Environment (MiAmbiente) of Panama within the framework of the ONUREDD program developed at the scale of 1:50,000 (Ministerio de Ambiente, 2012). Using the editing feature, we simplified the government generated LULC shapefile in ArcMap 10.8, merging several similar categories, and renamed the new groups to reflect the composition. The final categories based on the government LULC were mature broadleaf forest, secondary broadleaf forest, pasture, scrub/bushes, urban areas, banana plantations, mangroves, water bodies and sediment (Figure 2). Because more detailed land cover had been determined for the San San Pond Sak region, we manually edited the LULC shapefile, merging polygons and renaming them to reflect the vegetation map produced by Sjögersten et al. (2020) (Figure 2). The new “swampy” vegetation categories are flooded mixed forest, orey (Campnosperma panamense) dominated forest, raphia palm (Raphia taedigera) dominated, flooded vegetation (stunted forest and bog plain/saw grass), and mangroves (Phillips et al., 1997; Sjögersten et al., 2020). The “swampy” regions are located primarily within San San Pond Sak and along the coasts of Almirante Bay (Figure 2). The total area (km²) of each LULC type was calculated for the six focal catchments and the total Almirante Bay catchment (see Supplementary Table 3). LULC is unlikely to have changed significantly from 2013 to 2019 for the focal rivers and the Almirante Bay catchment, as minimal deforestation took place during this period (see Supplementary Table 4).

For statistical analysis, LULC categories were grouped into forest (mature and secondary forest), human impacted (pasture, scrub/bushes, urban areas, banana plantations), and swamp (flooded vegetation, mixed flooded forest, orey dominated forest, raphia palm dominated) (see Supplementary Table 5). Fluvial, which was not used in the land use statistical analysis, consisted of water bodies and sediment, and ranged from 0 to 1.5% of the total catchment area (see Supplementary Table 5). We used a partial least squares (PLS) analysis to determine which LULC types best explained the variation in each measured river constituent, since the river constituents are not independent (Carrascal et al., 2009). PLS analysis was implemented in JMP15, which uses minimization of root mean predicted sums of squares (PRESS) to identify the optimal number of factors (SAS Institute, 2019). Coefficients with a variable importance (VIP) >0.8 are considered important (SAS Institute, 2019). We examined linear regressions between each of the river water variables identified as important from the PLS analysis and percentage of each LULC type (forest, human influenced, or swamp). These models were used to predict the river physical and chemical characteristics for the entire Almirante Bay catchment based on LULC.

To understand how LULC has changed from 2000 to 2020, in ArcMap 10.8.1, we extracted the percentage of deforestation for the six focal catchments and the Almirante Bay catchment using a global dataset of forest deforestation (Hansen et al., 2013). We removed the pre-2000 banana plantation regions from our deforestation estimates as it was clear the changes were due solely to harvesting rather than deforestation.

The National River Monitoring Network

Panama's Autoridad Nacional del Ambiente created a baseline dataset of river water quality through the establishment of a National River Monitoring Network (NRMN) (ANAM, 2013). Two hundred seventy-six sites in 100 rivers within 36 major catchments were measured once during the wet season and once during the dry season from 2009 to 2012 (ANAM, 2013). Ideally, there were three NRMN sites on each river, one at high altitude, a second in the middle course of the river, and third near the river mouth or near the confluence in the case of a tributary.

Measured parameters included conductivity, temperature, pH, DO, total dissolved solids, turbidity, biological oxygen demand, fecal coliforms, nitrates, and phosphates. pH was measured in situ with an electrometer, conductivity with a conductivity meter, temperature with a thermistor, and turbidity with a nephelometer. Other parameters were analyzed in the ANAM (now MiAm) Water Quality Laboratory, including dissolved oxygen by acid modification, and nitrates and phosphates with a spectrophotometer (ANAM, 2013).

We assessed which NRMN sites were likely influenced by coastal processes. Any NRMN site with mean conductivity >1,013 μS/cm (salinity >0.5 ppt) was classified as under coastal influence if the following criteria were met for the site: (1) was located near the coast on flat, low-lying land, and (2) had at least one individual river sample per site with a conductivity >3,796 μS/cm (>2 ppt, maximum recommended salinity for agricultural use). Twelve percent of the NRMN sites were classified as coastal (see Supplementary Table 6).

A water quality index (WQI) developed by the National Sanitation Foundation USA (Brown et al., 1970; Dinius, 1987) and modified for Panamanian rivers (ANAM, 2003), ranged from 0 to 100 based on a weighted combination of nine measured physical and chemical parameters (ANAM, 2013). DO and fecal coliform were the critical WQI parameters and were given the largest weights (ANAM, 2013). A modified formula was used if any parameter besides DO and fecal coliform was missing, otherwise the sample was excluded (ANAM, 2013). Each NRMN site was placed into one of five categories based on mean WQI: highly contaminated (0–25), contaminated (26–50), slightly contaminated (51–70), acceptable (71–90), and not contaminated (91–100). Overall, 67% of the NRMN sites had acceptable water quality, 26% of the sites had slightly contaminated water, and 7.3% of the sites had contaminated water (see Supplementary Table 6).

To place the six focal rivers into a broader context, we used a PCA to compare them to the NRMN sites. This analysis used the physical and chemical parameters that both datasets had in common: conductivity, temperature, pH, and DO, ${NO}_{3}^{-}$ and ${PO}_{4}^{-}$ . Since nitrite is a small fraction of TON, ${NO}_{3}^{-}$ is comparable to TON (World Health Organization, 1978). For each river site we calculated the mean value for each variable for use in the PCA. One river site from the NRMN was not included in the PCA as it was an extreme outlier. The WQI for the NRMN was used to color code the points on a plot of PC1 vs. PC2 to show where in the multivariate space based on the subset of six variables, the sites with different quality fell.

Results

River Water Characteristics

The river water physical and chemical characteristics of the two swamp rivers were occasionally influenced by saltwater intrusions during the dry season (Feb-Mar; see below). The characteristics described below exclude this saltwater intrusion period. The median water temperature for each river ranged from 24.2 to 25.6°C (Table 1). The coldest (22.3°C) and warmest (28.5°C) reported temperatures were both recorded in Pastores (Table 1). The median specific conductivity for each river ranged from 58 μS/cm in Pastores to 565 μS/cm in the Canal with corresponding salinities of 0.03 ppt and 0.27 ppt, respectively (Table 1). The median pH for the rivers varied from 5.5 in Black to 7.6 in Changuinola (Table 1). The lowest pH recorded was 4.4 in Black and maximum pH was 8.4 in Pastores (Table 1; see Supplementary Figure 1C). The median values of DO for each river ranged from hypoxic at 2.2 mg/L in the swamp rivers, which never exceeded 4.8 mg/L, to a median of 9.2 mg/L in Changuinola (Table 1). Specific conductivity generally covaried temporally across the rivers, which included the seawater intrusion period (Table 2; see Supplementary Figure 1B) suggesting a regional driver. Temperature covaried temporally across all rivers, except when compared to Black, likely due to the seawater intrusion period (Table 2; see Supplementary Figure 1A). pH and DO values were temporally correlated only among some pairs of rivers supporting a role of more local or catchment-specific drivers (Table 2; see Supplementary Figures 1C,D). During the seawater intrusion at Black, the median specific conductivity was 44,095 μS/cm with corresponding salinity of 28.3 ppt (Table 1; see Supplementary Figure 1B). The lowest recorded DO concentrations occurred during the saltwater intrusions in Black where the median value was anoxic (0.05 mg/L) (Table 1; see Supplementary Figure 1D).

TABLE 1

Table 1. Summary results from river water physical and chemical constituent monitoring from the six focal rivers, with median, mean, standard deviation, minimum, maximum, and count.

TABLE 2

Table 2. Pairwise correlations (r) between the rivers across the fortnightly sampling of the six focal rivers^a.

The nutrient concentrations also varied between the rivers. The median TON concentrations had a large range, with the lowest values falling below the detection limit in Black to the highest value of 0.401 mg N/L in Nigua (Table 1). The median ammonium concentrations fell below the detection limit in Changuinola, Oeste, and Pastores and reached a maximum of 0.033 mg N/L in Black (Table 1). The median phosphate concentrations ranged from 0.014 mg P/L in the Canal to 0.096 mg P/L in Pastores (Table 1), with a maximum of 0.235 mg P/L in Oeste and a minimum of 0.004 mg P/L in the Canal (Table 1). TON, ammonium, and phosphate values were temporally correlated only among some pairs of rivers indicating a role of local or catchment-specific drivers (Table 2).

Analysis of Variance

Analysis of variance showed significant differences between many of the rivers for most of the variables measured. The strongest effects included the following. The two swamp rivers had significantly lower pH, DO, TON, and higher specific conductivity than the other rivers (Tukey post-hoc tests p < 0.05; Figures 4B–D, 5A). Nigua had significantly higher TON than the other rivers (Figure 5A). Ammonium levels in Black were significantly higher than levels in the Changuinola and Pastores (Tukey post-hoc test p < 0.05; Figure 5B), while none of the others differed statistically from each other. Pastores and Oeste had significantly higher concentrations of phosphates than the other four rivers (Figure 5C).

FIGURE 4

Figure 4. Box plots of river water constituents from the fortnightly dataset of the six focal rivers (see Supplementary Table 1) used in the statistical analyses for (A) temperature, (B) conductivity, (C) pH, and (D) dissolved oxygen. One-way ANOVA identified significant differences between the focal rivers (P < 0.001) (see Supplementary Table 8). Means with different letters are significantly different (P < 0.05) as per Tukey's post-hoc test. For comparison, data from the national river monitoring network (NRMN) in Panama (PAN) from 2009 to 2012 are included (ANAM, 2013). In the box plot the middle line is the median, and the box contains the interquartile range (50% of the data) and the whiskers extend to 25 and 75% of the data. X represents the mean.

FIGURE 5

Figure 5. Box plots of nutrient river water constituents from the fortnightly dataset (see Supplementary Table 1) used in the statistical analyses for the six focal rivers for (A) total oxidized nitrogen (TON = ${NO}_{3}^{-}$ + ${NO}_{2}^{-}$ , B) ammonium, and (C) phosphate. One-way ANOVA identified significant differences between the focal rivers (P < 0.001) (see Supplementary Table 8). Means with different letters are significantly different (P < 0.05) as per Tukey's post-hoc test. For comparison, the mean nutrient concentrations from national river monitoring network (NRMN) in Panama (PAN) rivers were measured from 2009 to 2012 (ANAM, 2013).

Multivariate Analyses of the Focal Rivers

Principal component analysis across our entire dataset shows that only 45.4% of the variation in the data can be described by the first two principal components (PC) (Figure 6A). DO, pH, TON and PO₄ increase along PC1 to the right, while temperature, salinity and NH₄ decrease in the same direction (Figure 6A). PC2 in contrast shows all factors other than DO, which contributes little to this axis, increase in the same direction. Visual inspection shows that the swamp rivers largely overlap each other, as do the mountainous rivers, but that the two groups can be largely distinguished along PC1 (Figure 6A).

FIGURE 6

Figure 6. (A) Principal component analysis (PCA) for all the river constituents. Principal components (PC) 1 and 2, the orthogonal linear combinations of standardized variables that maximize the variance of the observations (SAS Institute, 2021). Each observation is placed in the multivariate space defined by PC1 and PC2. (B) Quadratic discriminant analysis (DA) for all the river constituents. The axes (canonical1 & canonical2) summarize the between-category variation, to provide the maximum separation between the groups (i.e., each river). The multivariate mean for each group is shown by a plus sign.

Discriminant analysis after stepwise removal of non-significant factors retained DO, TON, ${PO}_{4}^{-}$ and pH. The results show that most of the rivers can be clearly distinguished (entropy R² = 0.89; −2LogLikelihood = 33.84; N = 86) (Figure 6B). The samples were attributed to the correct river 90% of the time. Black and the Canal were occasionally confused with each other (with 6/28 incorrectly classified as the other river), but were clearly separated from the four mountainous rivers (see Supplementary Table 7). Amongst the four mountainous rivers, Changuinola, and Nigua were classified correctly 100% of the time, but two Pastores samples were misclassified (with 2/17 incorrectly classified as Oeste; see Supplementary Table 7). This shows that the river water characteristics are distinct, as they can be distinguished statistically.

Land Use and Land Cover

LULC varied between the river catchments in our study. Changuinola had the highest forest land cover (90%). Pastores, Oeste, and Nigua had comparatively moderate forest cover and significant amounts of human-impacted land use (69:31, 64:36, and 54:45 forest: human-impacted, respectively; see Supplementary Table 5). The catchments of the Canal and Black (i.e., the swamp rivers), were almost entirely swamp (98 and 96%, respectively; see Supplementary Table 5).

The three main land cover types (forest, human impacted, and swamp) were good predictors of each of the river water characteristics in the partial least squares (PLS) analysis except for conductivity, for which a model with 0 factors was identified as minimization the root mean PRESS. This suggests that variation in conductivity/salinity measured in the swamp rivers is driven by factors other than land cover, most likely saltwater intrusion from the coast. For the remaining six characteristics, partial least squares models each resulted in the recovery of two land cover factors with VIP scores >0.8 (Figure 7), the threshold often used to determine which variables are important in explaining the data (Carrascal et al., 2009). Percent swamp cover was an important explanatory variable for all six characteristics, and percent forest cover was an important explanatory variable for four characteristics (temperature, pH, DO and ${NH}_{4}^{+}$ ) (Figure 7). The percent of human influenced land cover was an important explanatory variable for TON and ${PO}_{4}^{-}$ . PLS coefficients based on centered and scaled data show more swamp cover resulting in increases in temperature and ${NH}_{4}^{+}$ and decreases in pH, DO, TON and ${PO}_{4}^{-}$ (Figure 7). Higher forest cover acts in the opposite direction, with more forest cover resulting in decreases in temperature and ${NH}_{4}^{+}$ and increases in pH and DO (Figure 7). More human-influenced land use resulted in an increase in TON and ${PO}_{4}^{-}$ (Figure 7). The relationships between each water variable and each LULC type with VIP >0.8 are shown in Supplementary Figure 2. The relationship between each LULC type and the river constituent can also be illustrated in multivariate space using PCA analyses conducted on the average conditions in each river and one parameter of LULC (Figure 8).

FIGURE 7

Figure 7. Partial least squares (PLS) model coefficients for centered and scaled data for LULC (forest, human impacted, and swamp) for mean temperature, pH, DO, TON, ${NH}_{4}^{+}$ , and ${PO}_{4}^{-}$ , measured in the focal rivers (n = 6). Starred bars have a variable importance for projection (VIP) score >0.8, indicating importance of each land cover type for predicting each river constituent value (SAS Institute, 2021).

FIGURE 8

Figure 8. Principal component analysis (PCA) of the focal river means for seven water quality variables (A). The smaller panels include the seven water quality variables in the PCA, plus (B) percent human impacted cover, (C) percent forest cover, (D) percent deforestation in the last 20 years, and (E) percent swamp cover. The inset variable loadings are not to scale with the plot of points on the ordination axis, showing only relative magnitude and direction.

Although it was not possible to measure the river constituents for the contributing rivers of the entire Almirante Bay catchment (Figure 1), river water characteristics can be estimated based on the relationships between LULC and the mean values from of the six focal rivers (see Supplementary Figure 2). The LULC for the Almirante Bay catchment had 42.3% forest, 34.9% swamp, and 22.6% human influenced areas, and fell within the range of LULC encompassed by focal catchments (see Supplementary Table 5). Based on single variable regressions with LULC, the predicted mean physical and chemical characteristics for the entire Almirante Bay catchment are as follows: temperature at 25.2°C; pH at 6.8; DO at 5.4 mg/L; TON at 0.182 mg N/L; ammonium at 0.034 mg N/L; and phosphate at 0.059 mg P/L (see Supplementary Figure 2 and Supplementary Table 9).

Changes in LULC over the last 20 years have been low to moderate. From 2000 to 2020 there was 3.4% forest cover change (0.17%/year) in the entire Almirante Bay catchment (see Supplementary Table 4). Very little of this change occurred in the swamps, where there was <0.2% forest cover change in Black and the Canal, as they are within protected areas (Figure 2; see Supplementary Table 4). The bulk of this change occurred in the small mountainous catchments that drain directly into the Bay, with a decrease in forest cover by 9.0% in Nigua, 6.2% in Pastores, and 4.3% in Oeste (see Supplementary Table 4). Thus, in the past 20 years the physical and chemical characteristics of river water that contributes to Almirante Bay may have shifted slightly in response to these moderate land use changes.

Comparison With Nationally Monitored Rivers

Overall parameters for the focal rivers fell within the range of values reported for the NRMN dataset (Figures 4, 5). However, the median temperature for the focal rivers was 1.9°C cooler than the median temperature for NRMN sites (Table 1, Figure 4A; see Supplementary Table 6). There is the possibility that some of the differences between the two datasets (this study and the NRMN) stem from the two different monitoring protocols and from natural interannual variation. Black, which had saltwater intrusions, had similar conductivity to coastal NRMN sites (see Supplementary Figure 3B and Supplementary Table 6). The two focal swamp rivers had DO concentrations that were substantially lower than the coastal NRMN sites, but similar to contaminated NRMN sites (see Supplementary Figure 3D).

For the six focal rivers and the NRMN sites, the principal component 1(PC1) explained 49% of the variance (Figure 9). ${NO}_{3}^{-}$ and ${PO}_{4}^{-}$ were positively loaded and pH and DO were negatively loaded on PC1 (Figure 9). Generally, the acceptable water quality sites scored negatively, the slightly contaminated sites scored slightly positive, and the contaminated sites scored positively on PC1 (Figure 9). The PC2 explained 17% of the variance, where the positive loadings were related to high conductivity and warmer temperatures, including the coastal NRMN sites (Figure 9). There were three distinct NRMN site groups outside the central cluster of the PCA. They were as follows: (a) contaminated river sites which were almost universally located in Panama City with low pH, low DO and high nutrients (positive PC1 and negative PC2); (b) slightly contaminated coastal sites with high specific conductivity and high temperature (positive PC1 and positive PC2); and (c) acceptable sites located in the high mountains of Volcan Baru and Cerro Punto hillslopes with low temperature and high DO (negative PC1 and negative PC2; Figure 9). Four of the focal river sites were grouped in the central PCA cluster, with, Changuinola, Pastores, and Oeste, overlapping NRMN sites with acceptable water quality and the Nigua falling at the transition between acceptable and slightly contaminated water quality (Figure 9). The swamp rivers (Black and Canal) overlapped with category b above, consisting of slightly contaminated river sites with coastal influence (Figure 9).

FIGURE 9

Figure 9. Principal component analysis (PCA) of mean river constituents for the six focal river monitoring sites and for 275 sites within the national river monitoring network (NRMN). Principle components (PC) 1 and 2, the orthogonal linear combinations of standardized variables that maximize the variance of the observations (SAS Institute, 2021). Each observation is placed in the multivariate space defined by PC1 and PC2. Rivers from the NRMN are colored based on water quality index. River stations with coastal influence and the focal rivers are outlined in black. The value in brackets is the percentage of the variation accounted for by the PC. The inset with rays stemming from point 0,0 and indicates the weight (length) and direction of each variable on the two PC axes.

Discussion

River Differences and Similarities in Bocas del Toro

Although the focal river sampling sites are situated within 30 km of each other, the water quality of the rivers we surveyed differed significantly from each other. Salinity/conductivity covaried over time across the rivers suggesting that variation in these characteristics was driven by region-wide factors and are not specific to the different catchments, with temporal variation in part due to precipitation and evaporation. With the exception of Black due to periods of seawater intrusion, it is likely that temporal covariance in temperature between the rivers was due regional climatological variations in solar radiation and air temperature. In the swamp rivers, pH covaried overtime, with a large range in values, indicating that this region has lower buffering capacity and is sensitive to changes in biogeochemical processes and climate. The factors that were most instrumental in distinguishing the water samples from the six focal rivers from each other were pH, DO, TON and phosphate. These four factors covaried temporally only among subsets of the rivers suggesting more local drivers that differed between the different LULC groupings. Among the six focal rivers low pH, low DO, low TON separated the swamp rivers from the others. The relative concentration of ${PO}_{4}^{-}$ vs. TON effectively distinguished the mountainous rivers from each other.

Our multivariate analysis showed that four focal rivers have similar characteristics to NRMN sites, for which the WQI had been quantified. Our focal rivers had similar characteristics to NRMN sites with acceptable to slightly contaminated water quality. This is consistent with our observations of the local landscapes. Inferring from this relationship the focal rivers can be divided into four groups: swamp rivers with coastal influence and characteristics similar to rivers with slightly contaminated water quality (Black, and Canal); small mountainous rivers with characteristics similar to rivers with acceptable to slightly contaminated water quality (Nigua); small mountainous rivers with characteristics similar to rivers with acceptable water quality (Oeste and Pastores), and a large mountainous river with characteristics similar to rivers with acceptable water quality (Changuinola).

Coastal Peat Swamp Rivers

San San Pond Sak wetland, is a national protected area. Peat swamp rivers, like Black and the Canal, are unique ecosystems with distinct biogeochemical processes, yet they are understudied particularly in the Neotropics. One example of their unique characteristics is that the phosphate concentrations of the focal blackwater rivers are higher than the global average for natural rivers (Meybeck, 1982; Alkhatib et al., 2007). This phenomenon can be explained via the following mechanism. Anaerobic respiration causes denitrification, where river nitrates are reduced to N₂ gas which is then lost to the atmosphere (Knowles, 1982). When DO, nitrates, and nitrites are depleted, iron (III) oxide bound phosphate is reduced in sediments by dissolved sulfide releasing phosphate into the water (Chambers and Odum, 1990; Parsons et al., 2017). The high levels of phosphates in Black are particularly interesting because this river's catchment consists of a 8 m thick deep peat dome (Phillips et al., 1997) which is thought to receive nutrients only from precipitation, resulting in phosphate limitation in the terrestrial ecosystem (Sjögersten et al., 2011; Cheesman et al., 2012). More research is needed to assess the biogeochemical dynamics of the San San Pond Sak peat swamp forest and connectivity to the rivers.

Our results demonstrate that the standard WQI approach is not ideal for assessing the natural water quality of coastal rivers and peat swamp rivers. The loadings on PC2 in the analysis of the river sites within the NRMN show that many coastal sites had water characterized as “slightly contaminated,” with higher conductivity and higher temperature. While many coastal rivers have pollution runoff from upstream reducing the WQI, higher conductivity and higher temperatures are also natural characteristics of coastal rivers. The Canal and Black align with the coastal NRMN group, suggesting that they are slightly contaminated. The two swamp rivers in our study, however, are natural and undisturbed, with >95% swamp vegetation. They are typical blackwater rivers with low pH, low DO, and variable concentrations of nutrients (Vegas-Vilarrúbia et al., 1988; Küchler et al., 2000; Alkhatib et al., 2007; Horbe and Santos, 2009; Rixen et al., 2010; Ríos-Villamizar et al., 2014). The two focal swamp rivers likely align with “slightly contaminated” NRMN sites because they have DO concentrations approaching hypoxia. DO has the largest weight in the WQI calculation (ANAM, 2013). We did not measure fecal coliforms in the focal rivers, which is the other highly weighted variable. It is important, that care be taken in interpreting the WQI when it is applied to these peat swamp rivers and other similar coastal rivers. If the water quality of such rivers is characterized as slightly contaminated under the WQI but the catchment is primarily swamp or fecal coliforms are low, then reduced water quality is likely due to natural sources and should not be considered polluted. As saltwater intrusions can also give confusing results with respect to conductivity, until a WQI metric specifically designed for use with such coastal rivers is developed, interpretation of water quality indices in coastal swamp rivers should be made with caution.

LULC Impacts on River Water Characteristics

Our study focused on multiple catchments with varying LULC in the same geographical region in order to evaluate the relationship between catchment LULC and river water characteristics. We found results similar to published studies that have more frequently compared river chemistry before and after a change in LULC or upstream and downstream of a particular type of LULC (Brogna et al., 2017; Tanaka et al., 2021). For example, as the human influenced land use increased and forest cover decreased in the catchment, river nutrients increased and DO decreased. This was exemplified by the high TON in Nigua with a 54:45 forest:human-impacted land cover compared to the lower average TON in Oeste with 64:35 forest:human-impacted land cover. Previous work (Wickham et al., 2000; Wickham and Wade, 2002) has shown a non-linear effect of human impacted habitats on water quality, such that a small increase in urban areas can have a larger than expected overall impact on water quality. This may explain the seemingly disproportionately high levels of TON in the Nigua compared to the other mountainous rivers. This disproportionate effect may be due to the fact that urban areas, agricultural lands, and roads are often adjacent to waterways, as is the case in Nigua. A more extreme example is evident in the NRMN where river sites with the highest TON and ${PO}_{4}^{-}$ and lowest DO were polluted rivers in Panama City.

Our results allow the use of LULC to predict water quality for the entire Almirante Bay catchment both past and present. The small changes in LULC suggest that over the past 20 years the physical and chemical characteristics that contribute to Almirante Bay may have shifted slightly. Deforestation, although moderate in the Almirante Bay catchment, may have increased the TON concentrations as TON was closely aligned with forest cover. Reduced forest LULC may have slightly increased the overall river water temperature, ${PO}_{4}^{-}$ and ammonium levels, and slightly lowered the pH and DO. The peat swamp forests in San San Pond Sak are protected and little land cover change has occurred, thus it is likely that the contribution water to the Bay from the swamp is unchanged. The mountainous catchments, on the other hand, are likely to be important drivers of river water quality change in the Almirante Bay catchment. It has been suggested that land use change in the last 100 years has degraded estuarine water quality, as coral reef communities in Almirante Bay have declined (Cramer et al., 2012; Aronson et al., 2014).

Riverine Influence on the Almirante Bay Ecosystem

River water quality impacts estuaries around the global (Eyre, 1994; Nagy et al., 2002; Costanzo et al., 2003; Burford et al., 2012; Fernandes et al., 2018; Howley et al., 2018) with nutrient runoff contributing to coastal hypoxia and eutrophication (Wiseman et al., 1997; Gray et al., 2002; Howarth et al., 2011). Hypoxia is well-documented in Almirante Bay (Adelson et al., in review; Altieri et al., 2017; Lucey N. et al., 2020; Lucey N. M. et al., 2020; Weinstock and Collin, 2021). Our focal rivers, however, did not carry particularly high levels of dissolved nutrients when they passed our monitoring sites, and further investigation is needed to evaluate river nutrient fluxes into the Bay and to evaluate whether these contributions could influence hypoxia and/or ecosystem degradation.

Nutrient content in “pristine” rivers or streams is generally low, but they vary naturally depending on catchment characteristics, atmospheric sources, and geology (Stallard and Murphy, 2012). Therefore, it is difficult to assess the background levels of nutrients prior to human impacts and harmful anthropogenic effects (Mainstone and Parr, 2002). It is also tricky to compare to regulatory baselines, as such regional regulatory and quality standards for nutrients in rivers vary or are absent. However, natural levels of nitrogen in world rivers are around 0.101 mg N/L TON and 0.015 mg N/L ${NH}_{4}^{+}$ (Meybeck, 1982). The natural peat swamp rivers (Black and the Canal) have very low TON and double the ammonium compared to natural global rivers, but are typical of reducing and anaerobic environments (WHO, 2017). Nigua has around four times higher nitrate levels than those found in natural global rivers (Meybeck, 1982), potentially related to the amount of development in that catchment. However, overall, our rivers have much lower nitrate/nitrite concentrations than good/moderate (mean 2.5 mg N/L total) rivers in Europe (Poikane et al., 2019) and are much lower than drinking water limits at 11.9 mg N/L of TON (WHO, 2017).

Dissolved phosphorus is naturally present in rivers, with ${PO}_{4}^{-}$ around 0.01 mg P/L (Meybeck, 1982). Like nitrate, phosphate concentrations can span an order of magnitude at a single monitoring station (Meybeck and Moatar, 2012; Stallard and Murphy, 2012). All the focal rivers have higher phosphates compared to the average for natural global rivers (Meybeck, 1982), and two agricultural rivers in Puerto Rico (Stallard and Murphy, 2012). Several of our focal rivers have moderate to high forest cover, suggesting that geology, especially underlying the Oeste and Pastores, may be enriching the phosphorus levels in these small mountainous catchments. The phosphate levels in our rivers have similar or lower concentrations than good/moderate (mean 0.1 mg P/L) rivers in Europe (Poikane et al., 2019) and there is no international standard set for phosphates in drinking water (WHO, 2017).

Concentrations of dissolved nutrients in some of the focal rivers are elevated compared to idealized “pristine” conditions, and they would likely contribute to nutrient cycling in the Bay. However, assessing the potential impacts caused by river dissolved nutrient fluxes in the marine environment is challenging without additional information. First, river flow rates are necessary to calculate nutrient fluxes. Secondly, both dissolved organic forms and particulate forms (not measured here) contribute significantly to nutrient fluxes to the world's oceans (Meybeck, 1982) and are integral components of marine biogeochemical cycles (Seitzinger and Sanders, 1997; Berman and Bronk, 2003; Karl and Björkman, 2015). Finally, evaluation of terrestrial contributions solely through monitoring rivers overlooks groundwater and coastal town effluent runoff into the Bay (Mainstone and Parr, 2002; Valiela and Bowen, 2002; Bowen et al., 2007). With respect to this last point, the town of Almirante which is situated at the mouth of the Nigua, is likely a major source of anthropogenic nutrient addition to the back of the Bay. The proximity of the town of Almirante to the area of the Bay that experiences the most extreme hypoxia and temperature inversions (Adelson et al., in review) suggests a possible contributing role from town runoff, which is not captured in this survey of river nutrients. There are other towns however, in Almirante Bay, such as Bocas del Toro on Colon Island, where wastewater enters the Bay. The coast of these towns likely experiences more water exchange than Almirante due to their proximity to the entrance of the Bay and the open ocean (Adelson et al., in review). This suggests that anthropogenic sources of nutrients into the Bay, such as from the town of Almirante, may have a substantial impact on marine ecosystem health. Determining the nutrient fluxes into the Bay and identifying the relative importance of rivers (Rabalais, 2002), coastal towns (Silva et al., 2013), and the open ocean (Davis et al., 2014) will be a crucial step in the assessment and understanding of ecosystem degradation, eutrophication, and hypoxia in Almirante Bay.

Conclusions

Overall, we found that six focal rivers around Almirante Bay, in Bocas del Toro Panama have significantly different river water physical and chemical characteristics that can be linked to catchment LULC. Our results suggest that the peat swamp rivers had water quality that could be classified as slightly contaminated due to low pH and DO, but that these characteristics are naturally occurring in blackwater rivers rather than from pollution. Using the relationship between our river point measurements of water characteristics and catchment LULC, we developed regressions to estimate river constituent characteristics for the entire Almirante Bay catchment. Deforestation rates although low to moderate in the Almirante Bay catchment in the last 20 years, may influence river water quality by increasing TON inputs from the mountainous catchments. This study shows that non-point source nutrient input, from wastewater and urban runoff from coastal towns, may have a greater impact on coastal hypoxia than river runoff.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

KC and RC designed the study with input from SG, GP, KD, and MT. VB, XB, CC-Á, and KC collected the samples and data in the field with occasional help from AA. KC managed the data QA/AC, conducted the GIS analysis, and generated the maps. KC and RC analyzed the data and wrote the manuscript. SG, GP, KD, MT, and AA commented on the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the following grants NSF: (OCE - 1924220 to RC), (OCE - 1924551 to GP, SG, and AA), and (OCE - 1924664 to KD). KC and AA were partially supported by fellowships from STRI.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We are grateful to the Ministerio de Ambiente de Panama (MiAmbiente) for permission to work in Panama and for providing access to the national river monitoring network results. We thank T. Romero, A. Bielnicka, and the staff of the STRI soils laboratory for conducting the nutrient analysis, P. Góndola and the staff of the Bocas del Toro Research Station for logistical support, F. Bellido, E. Brown, D. González, and H. Arlauskas for assistance in the field, A. Ibáñez and E. Arauz for editing the Spanish-language abstract.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frwa.2022.737920/full#supplementary-material

References

Alkhatib, M., Jennerjahn, T. C., and Samiaji, J. (2007). Biogeochemistry of the Dumai River estuary, Sumatra, Indonesia, a tropical black-water river. Limnol. Oceanogr. 52, 2410–2417. doi: 10.4319/lo.2007.52.6.2410