Influences on tidal channel and aquaculture shrimp pond water chemical composition in Southwest Bangladesh

Detailed geochemical studies of both major and minor elements in Bangladesh surface waters are sparse, particularly in shrimp aquaculture pond environments. Therefore, water samples from shrimp aquaculture ponds and tidal channels were collected in high precipitation (July) and low precipitation (May) months from 2018–2019 in Southwest Bangladesh and analyzed for complete water chemistry. Selenium (Se) and arsenic (As) were elevated above WHO guidelines in 50% and ~ 87% of samples, respectively, but do not show any recognizable spatial patterns. Shrimp pond and tidal channel water compositions in the dry season (May) are similar, illustrating their connectivity and minimal endogenous effects within shrimp ponds. Tidal channels are less saline in July than shrimp ponds still irrigated by tidal channels, suggesting that either farmers limit irrigation to continue farming saltwater shrimp, or the irrigation flux is low and leads to a lag in aquaculture-tidal channel compositional homogenization. δ18O and δ2H isotopic compositions from samples in May of 2019 reveal tidal channel samples are closer to the local meteoric water line (LMWL) than shrimp pond samples, because of less evaporation. However, evaporation in May shrimp ponds has a minimal effect on water composition, likely because of regular drainage/exchange of pond waters. Dissolved organic carbon (DOC) is positively correlated with both δ18O and δ2H in shrimp ponds, suggesting that as evaporation increases, DOC becomes enriched. Multiple linear regression reveals that As and Se can be moderately predicted (adjusted R2 values between 0.4 and 0.7, p < 0.01) in surface waters of our study with only 3–4 independent predictor variables (e.g., Ni, V and DOC for Se prediction; Cu, V, Ni and P for As prediction). Thus, this general approach should be followed in other regions throughout the world when measurements for certain hazardous trace elements such as Se and As may be lacking in several samples from a dataset. Supplementary Information The online version contains supplementary material available at 10.1186/s12932-021-00074-2.


Introduction
Although there has been much research in Bangladesh on groundwater and contaminants such as arsenic (As) (e.g., [4,7,17,38,47,48]), less emphasis has been placed on surface water chemistry, especially in Southwest Bangladesh. Multiple studies in Bangladesh have geochemically examined tidally dominated rivers (hereafter called tidal channels) and adjacent waterways in both temporal and spatial ways (e.g., [2,3,24,33,59,63]). Studies have also looked at trace element concentrations in aquaculture ponds near the coastal region of Bangladesh [23,32,62] or in coastal Bangladesh rivers [36,52] and nonspecified lakes/ponds [1]. However, these studies have focused more on risk assessment, general reporting of trace element or major element concentrations, and overall water quality. While reporting concentrations and changes in water chemistry is useful for preliminary work, more detailed understanding of surface water geochemical relationships is imperative for predicting where waters may have elevated concentrations of certain hazardous elements such as As or selenium (Se), and what ultimately controls the element concentrations.
Studies thus far also have not thoroughly researched the relationship of tidal channel waters with aquaculture shrimp ponds, which are often irrigated with tidal channel water during the dry season in Southwest Bangladesh instead of groundwater [6] (Fig. 1). Recent studies indicate that tidal channel waters often have elevated arsenic, especially in the dry season [6] [25], so there is concern that shrimp ponds and the shrimp grown in them may also have high arsenic.
Particularly reliant upon surface water in Southwest Bangladesh (e.g., Satkhira and Khulna districts) are aquaculture and agriculture, where rice farming, fishing, and shrimp farming are the primary sources of income for people in the area (e.g. [11]) and shrimp farming in particular has increased dramatically since the 1980s [9]. Shrimp ponds (bottom left of Fig. 1) are typically 1 m or less in depth on average and vary largely in area, although averaging around 3 ha, while the tidal channels irrigating them have variable depths ranging up to several meters.
In general, critically understanding relationships between trace elements in surface waters in Bangladesh is particularly needed because trace element chemistry in Asian rivers is poorly understood [28], and the aquaculture shrimp ponds and tidal channels in Southwest Bangladesh are an innate offshoot of the large Ganges River system. Additionally, better understanding of water chemistry in Bangladesh aquaculture environments has worldwide applicability to coastal aquaculture systems in other countries such as India and Taiwan. Thus, this study aims to holistically examine possible influences (i.e., precipitation, evaporation, irrigation) on the composition of tidal channel and shrimp pond water in Southwest Bangladesh, and to use multivariate techniques to assess whether As and Se concentrations can be predicted in surface waters. The main research questions are: (1) Do early monsoonal rains affect water compositions of shrimp ponds and tidal channels differently? (2) Are any elements or other geochemical parameters sensitive to endogenous changes within shrimp ponds following tidal channel irrigation? and (3) Can Se and As surface water concentrations be reasonably explained and modeled by a small subset of other variables through multiple linear regression modeling?

Study area
Southwest Bangladesh is a tidally influenced area near the coast of the Bay of Bengal (Fig. 2), with the tidal influence extending just north of Khulna City in the dry season (e.g., [16]). Studies in the area have shown that tidal influence on tidal channel water composition varies by season and the amount of freshwater discharge from Fig. 1 Conceptual diagram of the water cycle within Southwest Bangladesh and example images of tidal channels that irrigate shrimp ponds (both a typical main channel like that sampled in this study (A) and an inlet from a smaller, connective channel (B). A typical shrimp pond (C). A video of tidal channel irrigation to a shrimp pond is provided as Additional file 2 upstream [53,54]. Specifically, salinity ranges between 0 and 2 parts per thousand (ppt) between the tides, and pH ranges 0-0.4 units between tides throughout the year depending on location within the Sundarbans [53,54]. This tidal area has not been as extensively studied as the northern floodplains in the G-B-M delta or the Himalayan foothills, where river water and sediment compositions were previously researched (e.g., [14,43,61,64]). The tidal region includes part of the Sundarbans mangrove forest, although much of this forest was converted to agricultural islands back in the 1960s-1970s (e.g., [51] and the references cited therein). These agricultural islands are also known as "polders, " and are surrounded by embankments, shielding these islands from storm surges or tidal inundation [5,51].
Bangladesh experiences a strong monsoonal climate, where biseasonal precipitation causes approximately 80% of yearly rainfall to occur between the months of June to September (e.g., [11,18,51]). This large seasonal difference in precipitation can lead to vast changes in surface water chemistry during different seasons, particularly in coastal regions (e.g., [2,59,63]).

Sample collection
Surface water samples were collected in both the early monsoon season in July of 2018 and 2019, and during the end of the dry season in May of 2019. July average precipitation in Khulna is ~ 363 mm, while May average precipitation is ~ 180 mm [10]. All samples were taken within the general vicinity of Southwest Bangladesh, and predominantly within the Khulna District (Fig. 2). All but one shrimp pond aquaculture sample (KA-10) were irrigated with tidal channel water.
Samples were collected as follows: • July: 11 shrimp ponds + 5 tidal channel samples + 1 rainwater sample (n = 17) • May: 10 shrimp ponds + 4 tidal channel samples (n = 14) May and July shrimp pond samples were taken in the same relative locations to avoid differences in tidal influence and thus allow more direct comparison. Tidal channel samples were taken at different locations between seasons, but those used for dissolved ion concentration comparison in spider diagrams (P32-TC, P32-1 and P32-2 in July; MD-TC-16, MD-TC-18, MD-TC-19, and MD-TC-22 in May) were collected in channels that irrigated the studied shrimp ponds in both seasons. One rainwater sample was collected in July of 2018 from a tin roof collection device.
Surface water samples were collected with plastic syringes and buckets from the upper 1 m of the water column and were rinsed at least once between each sample. Samples were filtered in the field with a syringe and 0.45 µm polypropylene filter into 60-and 125-mL bottles. May 2019 and July 2019 samples were all filtered into 60 mL bottles until the bottles were completely filled to prevent oxygen and hydrogen isotope exchange with headspace before analysis. Only 2019 samples were collected for isotopic analysis because in our lab, analytical capabilities for oxygen and hydrogen isotopes were not yet available in 2018. The July 2018 tidal channel and rainwater samples were placed in 500 mL bottles, with the tidal channel sample filtered in the laboratory with a 0.45 μm polypropylene filter (not filtered in the field).

Field measurements
A Hach Hydrolab DS5 was used to gather in situ surface water measurements of pH, oxidation-reduction potential (ORP) in millivolts (mV), specific conductivity (SpC) in microsiemens per centimeter (µS/cm), and temperature in degrees Celsius (˚C) for July 2018 shrimp pond, May 2019 shrimp pond, and May 2019 tidal channel samples following the methods of Ayers et al. [6]. July   2018 tidal channel SpC was measured with a hand-held  SpC Hanna probe while one July 2018 shrimp pond sample (KA-1) SpC value was estimated through total ion  summation. ORP measurements from the Hach Hydrolab DS5 were used for distinguishing oxic and anoxic conditions for surface water samples, as described in [6,7]. ORP (relative to the Ag/AgCl redox couple) was converted to Eh (relative to the standard hydrogen electrode (SHE)) in Table 1 by adding 187 mV to each field measurement.

Sample preparation and analyses
2-10 mL of each filtered sample was acidified with a drop of concentrated nitric acid (HNO 3 ) in the laboratory for inductively coupled plasma (ICP) analysis, while approximately 20-30 mL of sample was left unacidified and used for ion chromatography (IC) and total organic carbon analysis (TOC). Approximately 2 mL of filtered, unacidified May and July 2019 samples were transferred to glass vials via plastic pipette immediately upon opening the 60 mL bottles for δ 18 O and δ 2 H isotope analysis.
All acidified water samples were analyzed on a Perkin Elmer NexION 2000B ICP-MS in both standard and kinetic energy discrimination (KED) modes using EPA Method 6020B at Vanderbilt University for the elements As, Be, Cd, Co, Cr, Cu, Fe, Mn, Mo, Ni, Pb, Sb, Se, Ti, Tl, V, and Zn. All acidified surface water samples were also analyzed on an Agilent 5110 VDV ICP-OES using EPA Method 6010D at Vanderbilt University to report the ions: Al, As, B, Ba, Ca, Fe, K, Mg, Na, P, S, Si, and Sr. All 2018 filtered, unacidified surface water samples were analyzed for inorganic and organic carbon content with a Shimadzu model TOC-V CPH/CPN using ASTM Method D-7573-09 at Vanderbilt University, as described in [6,7] USGS50 were used to correct measured sample values through the USGS LIMS for Lasers data reduction scheme.
More details on sample preparation are included in the Additional file 1.

Quality control
Method blanks of samples were taken with deionized water in the field in July 2018 and May 2019 and analyzed for complete chemical analysis (IC, TOC, ICP-MS, ICP-OES), and routinely show concentrations at or below method detection limits ( Table 2). For samples run on the ICP-MS, ICP-OES, IC and TOC analyzers, periodically measured concentrations in standards were required to be within 15% of the known value and blanks were required to be below the method detection limit (MDL). A duplicate shrimp pond sample from May 2019 showed all concentrations except Mn in very close alignment with one another (3.0% average difference in concentrations for all reported elements (including DIC and DOC) that had concentrations > 10 ppb), illustrating the representative nature of each sample. Method detection limits are listed in Table 2. Average charge-balance error was 4.6%, similar to the average charge-balance error in [6] of 3.9%.

Data reduction
Maps were generated using ArcGIS 10.6.1 and QGIS 3.10, while RStudio and Microsoft Excel were used for figure generation and statistical analyses. Multiple linear regression subset selection used the "leaps" package in R, with 8 variables used as the max number of subsets [42]. The Geochemist's Workbench 14.0 was used to generate a Piper diagram and an evaporation model in the React program, using sample MD-TC-22 as the input basis. HCO 3 − concentrations (from measured dissolved inorganic carbon (DIC)), charge imbalance, and mineral saturation indices were calculated using SpecE8 software in The Geochemist's Workbench 14.0 with the default thermo.dat database [13]. Uncertainties are reported as sample standard deviation (1σ).

Geochemical concentrations
Dissolved concentrations tend to be lognormally distributed, so they are summarized using the geometric mean. Element concentrations reflect relatively saline waters in shrimp ponds regardless of season, shown by high conservative ion concentrations approaching average seawater concentrations ( Table 2). Tidal channel waters show more variability in conservative ion concentrations between sampling months compared to shrimp ponds. Eh values in both sampling months are generally consistent and positive in surface waters, indicative of oxidizing conditions (422 ± 48 mV). Including both sampling months (because of minimal seasonal variability in DOC), the mean shrimp pond DOC concentration is 7.3 ± 4.3 mg.L −1 , while the mean tidal channel DOC concentration is 3.1 ± 1.0 mg. L −1 , compared to a river world average of 5.8 mg. L −1 and 12 mg. L −1 median value in the world's eutrophic lakes [6] and the references therein). pH values show relatively little variability regardless of sampling month, although shrimp pond waters at 8.4 ± 0.4 are markedly more alkaline than tidal channel waters at 7.6 ± 0.1. Out of elements with known adverse health outcomes in excessive quantities, Cr (6.6 ± 3.6 μg. L −1 ) and Mn (3 ± 17 μg. L −1 ) overall surface water values are all well below levels of health concern [74], while Pb, Be, Cd, and Tl were either close to detection limit or below detection limit and thus not reported. Ti, Sb, Fe and Mo values are also not reported because of either low values (e.g., nearly all May samples had negative reported Fe values) or concerns with interference on the ICP. However, both As (25 ± 19 μg. L −1 ) and Se (45 ± 45 μg. L −1 ) are at levels of potential health concern across both sampling months based on WHO and EPA drinking water guidelines, which are 10 μg. L −1 for As for both the EPA and WHO, and 50 μg. L −1 and 40 μg. L −1 , respectively, for Se EPA and WHO guidelines [70,74].

Element and geochemical parameter correlations
Concentrations of most conservative elements are strongly correlated, indicative of their similar mobilities in solution (Additional file 1: Fig. A1). Concentrations of many elements that tend to behave nonconservatively     All measurements in mg.L

Monthly composition differences and element enrichment relative to seawater
In general, lower element concentrations are seen in July tidal channel waters relative to May tidal channel waters (Fig. 3). Concentrations in July tidal channel waters, however, show much more variance, although nonconservative elements have less consistency in variation. July and May shrimp ponds show slight differences in average compositions, albeit with most overlap contained within 1σ error bars (Fig. 4a). However, Se, Cu, Cr, and Co are all elevated in July shrimp ponds relative to May shrimp ponds outside of 1σ variation. July shrimp ponds are slightly less saline on average than May shrimp ponds (Additional file 1: Fig. A3). When normalizing shrimp pond (Fig. 4b) and tidal channel (Additional file 1: Fig.  A4) concentrations to that of seawater (values from literature; Additional file 1: Table A1), several elements are clearly enriched, such as Ba, As, Se, and Zn. Additionally, DOC is well enriched in shrimp ponds relative to values typically seen in the Indian Ocean [50], shown through DOC plotting well above 0 in Fig. 4a.

Tidal channel versus shrimp pond composition
May tidal channels and shrimp ponds have similar average compositions, although there is more variability in the nonconservative elements such as Zn and Ni (Fig. 5a). In July however, tidal channel and shrimp pond compositions are not as similar, although it is noted that there is high variability within July tidal channel samples (Fig. 5a). In general, July shrimp ponds have higher average concentrations of elements and DOC compared to July tidal channels.
Oxygen and hydrogen isotopes specific trends δ 18 O and δ 2 H from tidal channel samples (Table 3) plot relatively close to the local meteoric water line (LMWL) (Fig. 6). Furthermore, they plot linearly between the LMWL and the expected value for seawater. Shrimp ponds plot farther off the LMWL than tidal channel samples and their stable isotope values are heavier (more positive) than tidal channel isotope values. δ 18 O and δ 2 H show a strong positive correlation with DOC in shrimp ponds (δ 18 O and DOC; Spearman ρ = 0.92, p = 0.00047) (Fig. 7). As shrimp ponds become more enriched in both 18   The model estimates that only a relatively small fraction of water evaporation (~ 0.5-10%) is needed to concentrate non-reactive major conservative ions in solution to the concentrations observed in the irrigated shrimp pond (KA-15). When looking at modeled salinity (essentially grouping all major seawater ions together), 1.5-2% evaporation of water can explain the relative increase in salinity observed from tidal channel to shrimp pond  Fig. A5), which is close to intersample variability.
In addition to the aforementioned mass-balance focused geochemical model, relatively straightforward isotopic equations can be used to estimate the amount of evaporation in a shrimp pond after tidal channel input. Using several simplifying assumptions such as thermodynamic equilibrium and an equilibrium fractionation factor (α (l/v)) of 1.0098 at 20˚C (even though sampling temperatures were ~ 30-35˚C) for 18 O/ 16 O exchange [44], a Rayleigh distillation equation may be used to estimate the residual fraction of water (f ) left in a system (Eq. 1), [29]: where R is the isotopic ratio of the residual water and R 0 is the isotopic ratio of the initial water. In a system where water is evaporating and the liquid phase is known to be more enriched in the heavier stable isotope (i.e., 18 term for the fractionation factor (α) becomes < 1 (1/α), causing isotopic enrichment in the residual liquid phase as f decreases. After converting the δ 18  . However, this also assumes no other influences on stable isotopes in the shrimp pond, such as precipitation, and assumes that all water was derived from the tidal channel. These are reasonable assumptions, as during the sampling period in May 2019 no measurable precipitation was observed, and the impermeable clays lining most of the ponds prevents groundwater exfiltration or seepage from other water sources [6,7].
The discrepancy between the estimated 2% evaporation based on salinity and 10% evaporation based on isotopes may be in part due to slight non-conservative processes occurring in ponds such as ion exchange that affect saltwater ions, and our isotopic evaporation model is merely . Less similarity is seen in July (B), with large variance seen in July tidal channels and lack of overlap between the two water types. Two Si values are omitted from the July shrimp pond data because of anomalous negative values an estimate, as we assume an open system with Rayleigh fractionation. In reality, the system is likely somewhere between an ideal open and ideal closed system. Additionally, fractionation is likely slightly different at the warmer temperatures (> 20˚C) that occurred during our study.

Selenium and arsenic spatial and seasonal variability
Arsenic concentrations in shrimp ponds and tidal channels are close to or above WHO drinking water guidelines (10 μg. L −1 ) in both May and July samples throughout SW Bangladesh (Additional file 1: Fig. A5). There does not appear to be any spatial correlation between measured arsenic values in both shrimp ponds and tidal channels. Additionally, selenium values are close to or above WHO drinking water guidelines (40 μg. L −1 ) throughout SW Bangladesh in tidal channels and shrimp ponds, with seemingly random spatial concentration trends (Additional file 1: Fig. A7). In general, spatial heterogeneity in  surface water trace element compositions appears to be the norm. While both Se and As are within 1σ variation in May and July tidal channels (Fig. 3), Se is greater in July shrimp ponds compared to May shrimp ponds outside of 1σ variation (Fig. 4a). Arsenic is within 1σ variation when comparing May and July shrimp ponds (Fig. 4a).

Multiple linear regression
Best subsets multiple linear regression (details on variable selection in the Supplementary Document) suggests that the variables Cu, P, V and Ni can explain most As variance (Fig. 8), with an adjusted R 2 of 0.39 (p = 0.006) (Additional file 1: Table A3). Their significance is supported by V, Ni, and Cu each having p-values < 0.05 when looking at those four variables together as a linear regression fit for As (Additional file 1: Table A3). For Se, linear regression using the concentrations of only V, Ni, DOC, and Cl result in one of the best adjusted R 2 values (0.74) (Additional file 1: Table A4; Fig. 8), and each variable has a p-value < 0.01 (Additional file 1: Table A4). When removing Cl from the predictive model, the adjusted R 2 becomes 0.67 (p < 0.0001) with each variable (V, Ni and DOC) having a p-value < 0.01 (Additional file 1: Table A5).
Details of other models and violations of multiple linear regression assumptions are provided in the Additional file 1.

Trace element enrichment in surface waters
Trace elements of initial focus in this study were Mn and As, because of their known elevated concentrations in groundwater in the region (and potential to enter surface water through groundwater irrigation upstream or exfiltration) and known adverse health effects in excessive quantities (e.g., [27,34,56,58]). Furthermore, [6] measured elevated As > 10 μg. L −1 (WHO drinking water guideline, [74]) in 78% of shrimp ponds on Polder 32 in SW Bangladesh (several samples in this present study are from/around the same polder), with 71% of May tidal channels around Polder 32 containing As > 10 μg. L −1 as well. [6] also measured Mn > 400 μg. L −1 (WHO healthbased value, [74]) in 6% of surface water samples.
In this study, however, no dissolved Mn concentrations were above the WHO health-based value of 400 μg. L −1 [74] in either shrimp ponds or tidal channels, which is expected based on Mn solubility decreasing when pH is > 6, due to the formation of solid-phase Mn oxides and hydroxides (e.g., [35]). Thus, elevated Mn values were reported in some of [6]'s shrimp pond samples despite similar pH and redox conditions. Arsenic was elevated above the WHO standard in almost all shrimp pond and tidal channel samples regardless of season, with ~ 87% containing As > 10 μg. L −1 ( Table 2). This supports data from [6] where elevated As in SW Bangladesh surface waters was also reported. However, dissolved As is not expected to be as elevated in near-neutral pH, oxidizing conditions such as those seen in surface waters (e.g., [6,65]), and thus more research on As cycling at the groundwater-surface water interface in Bangladesh is needed (e.g., [12]).
Arsenic values were > 10 μg. L −1 in 75% of July samples (shrimp pond and tidal channels) during the early monsoon in this study, while [6] measured only 11% of tidal channels in October with As > 10 μg. L −1 . This suggests that the monsoon rainfall has not yet reached a high enough continuous level in July to thoroughly dilute the source of surface water As compared to October. In Khulna the average precipitation in the preceding four months were ~ 1,261 mm in October, only ~ 142 mm in May, and ~ 627 mm in July [10].
Because there were no previous reports of high Se in SW Bangladesh surface waters, it was not an initial element of concern going into this study. Furthermore, Se has been reported as low in Bangladesh groundwater Fig. 8 Multiple linear regression models for grouped variables initially selected for Se and As, with adjusted R 2 values along the y-axis. The best model was computed for each subset size (e.g., each amount of predictor variables utilized in the overall model). n = 25 because of some samples missing predictive variable measurements and thus being omitted, which also explains the slight discrepancy between multiple linear regression analyses completed with more samples (i.e., Additional file 1: Table A4) (e.g., [26,27,55]), so even if heavy groundwater irrigation or exfiltration could affect surface water concentrations, Se was not expected to be elevated. However, dissolved Se was found to be > 50 μg. L −1 (EPA MCL; [70]) in 50% of all samples and > 40 μg. L −1 (WHO guideline; [74]) in 50% of all samples (Table 2), with substantially higher average concentrations in the July shrimp ponds compared to May shrimp ponds (Fig. 4a). While low levels of Se intake can lead to nutritional deficiencies in humans such as impaired growth, or thyroid function abnormalities (e.g., [73]), excessive Se intake can be toxic for humans and animals, particularly for inorganic species of Se [68]. However, when shrimp ponds are often converted to rice paddies in the wet season in SW Bangladesh, the chemically reducing paddy sediment likely leads to the formation of insoluble elemental Se [73].
When normalized to average concentrations in seawater, As, Se and Mn are all enriched in shrimp ponds (Fig. 4b), while As and Se are also enriched in tidal channels (Additional file 1: Fig. A4). Thus, the sourcing of these "enriched" elements cannot be explained by seawater from the Bay of Bengal alone, and must have outside sources such as riverine flow, direct anthropogenic input (fertilizers or other chemical supplements), or groundwater flow. The higher Se concentrations in July shrimp ponds versus May shrimp ponds is difficult to explain (Fig. 4a), particularly because the same distinct seasonality is not seen in tidal channels (Fig. 3). A possible explanation is that desorption mechanisms are at work in the July shrimp ponds, although there are no correlations between Se and Eh or pH to suggest this is redox or pH influenced (Additional file 1: Fig. A2). Because colloids can have an important influence on trace elements in solution (e.g., [28]), future spectroscopy and speciation analyses may help identify whether they play an important role in enrichment and seasonality in Bangladesh.

Tidal channel and shrimp pond connectivity
Strong compositional similarities were seen between May shrimp ponds and the tidal channel sources irrigating them (Fig. 5a). This suggests that there is little compositional change within the shrimp ponds after irrigation, although there is more variance between nonconservative elements, suggesting that sorption/desorption reactions may be occurring after shrimp ponds have been irrigated. This may particularly be true for cation species like Ni and Zn, which may undergo desorption in more saline waters due to competition for sorption sites with salt cations such as Ca and Na (e.g., [60,66]). Across all dissolved water samples in this current study, moderate positive correlations between metal cations such as Co, Ni, and Cu (Additional file 1: Fig. A1) suggest similar geochemical processes occurring in the natural environment, such as sorption/desorption effects. July shrimp ponds show much less compositional overlap with July tidal channels (Fig. 5b). However, there is high variance of dissolved ion concentrations within July tidal channel samples (Fig. 3), which could be indicative of salinity fluctuation closely tied to rainfall and river discharge. July shrimp ponds and tidal channels may not be as similar compositionally because of shrimp farmers limiting irrigation of less saline waters to prolong harvest of brackish water shrimp varieties, or because of slow mixing rates between the irrigation source (tidal channels) and the shrimp ponds, which may not become complete until late in the monsoon season (i.e., late August or September). Future work analyzing water fluxes in and out of shrimp ponds over time would help better quantify the mixing of pond water and irrigation source.

Enrichment of As and Se
Arsenic and selenium concentrations vary widely throughout the region. However, they are both ubiquitously elevated > WHO drinking water guidelines in May and July shrimp ponds and tidal channels (Additional file 1: Figs. A6, A7). The widespread occurrence of high As and Se concentrations throughout the region suggests that As and Se are predominantly not locally point-sourced.
Arsenic may be sourced from arsenic-rich groundwater in the region (often > 10 μg. L −1 ) as previously suggested [6], upstream weathering in the Himalayas, or anthropogenic effluent run-off from cities upstream such as Khulna. Because As shows little to no correlation with conservative ions found in seawater and tends to behave nonconservatively (Additional file 1: Fig. A1), it is difficult to deduce whether As increases as more groundwater is present in the system, particularly because there are two possible lower salinity end-members (riverine flow from upstream and groundwater flow) with potentially different As concentrations. It is noted that a January (dry season) freshwater river sample in the nearby Meghna River (where groundwater As concentrations are > 10 μg. L −1 ) had As < 10 μg. L −1 [12], although this river has a different drainage basin.
Arsenic concentrations in this study are also much greater than surface waters sampled in 2016 in the same geographic region [8]. The differences in As concentrations between [8] and both this study and [6] may be attributed to small sample size, overall less saline (and thus potentially more diluted) samples in Ayers et al. [8], and sample location heterogeneity, as all three studies used similar methodology. For example, if groundwater exfiltration is a possible source of As to surface waters, sampling locations from [8] may receive different exfiltration rates of groundwater because of widespread geologic heterogeneity in the area (e.g., [7]).
Selenium concentrations higher than average seawater in estuaries can be caused by conservative mixing of polluted/enriched Se river water and relatively Se-depleted seawater, such as in the Solent area in mid-southern England [46] or the Rhone river delta off the Mediterranean Sea [31]. Anthropogenic input of Se was documented in the San Francisco Bay estuary from both sewage treatment plants and oil refineries, albeit the concentrations in the mixed estuary samples were much lower than those seen here [21,22]. Selenium was also elevated relative to natural background in agricultural irrigation drainage in central California in the mid-1980s [40]. Thus, Se may be sourced from sewage, oil refineries, or agricultural runoff in the area. It is also interesting that Se concentrations are much higher in our samples from 2018 and 2019 (regardless of season) compared to results from surface waters in the same general study area in 2016 analyzed using the same methods [8]. One possible explanation for this is that Se is anthropogenic and the source of its release into the dissolved load began after sampling in 2016, potentially from industrial centers in Khulna [25], while another possible explanation is that the overall less saline samples taken in 2016 were more diluted with rainwater or other sources than ours in 2018 and 2019, thus lowering Se concentrations.
Recently, in the Ganges River system, "hot spots" of elevated trace element concentrations (relative to "background" in the study) were observed near large urban and industrial areas, but these concentrations became diluted by other river tributaries downstream [15]. Thus, a city like Khulna could be a source of Se or other trace elements, and dilution is limited until either the tidal channels empty into the Bay of Bengal or it becomes peak monsoon season. Regardless of the source, a dramatic increase in Se input would be necessary for an increase in surface water concentrations from ~ 0.5-2 μg. L −1 in 2016 to ~ 10-200 μg. L −1 in 2018-2019.
Surface water isotopic composition, trends, and evaporative effects δ 18 O and δ 2 H isotopes in surface waters can be useful proxies for illustrating mixing or evaporative processes. Any samples plotting on the local meteoric water line (LMWL) are theoretically local precipitation, while samples plotting to the right of the LMWL have undergone fractionation, which can be interpreted as evaporation [29], or mixing with isotopically enriched water such as seawater. All tidal channel samples plot close to the LMWL and trend towards the isotopic value of seawater as they become heavier or more enriched, suggesting they have undergone minimal evaporation and represent mixing between rainwater and seawater (Fig. 6). Dry season (May) shrimp pond samples are heavier in δ 18 O and δ 2 H and plot farther off the LMWL compared to tidal channel samples, suggesting they have undergone more evaporation and are not simply depicting mixing between rainwater and seawater.
However, even though shrimp ponds are likely undergoing more evaporation, there is not a positive correlation with isotope values and salinity (SpC) (Additional file 1: Fig. A8). A positive correlation would be expected, because as water evaporates, conservative ions should stay in solution and become more concentrated. The lack of an observed positive correlation is likely a consequence of each shrimp pond having different starting salinities following irrigation from tidal channels, with concentration variations in the tidal channel irrigation source "masking" the relatively minor effects of evaporation on conservative salt ions. However, when looking at a shrimp pond and irrigating tidal channel directly adjacent to it, slight evaporative enrichment occurs with conservative elements, while several nonconservative elements are depleted in shrimp ponds (i.e., Mn and P; those less likely to desorb from saltwater cation exchange), possibly due to mineral precipitation or sorption onto sediment surfaces (Additional file 1: Fig. A9).
Evaporative models further illustrate the aforementioned point that although evaporation is likely occurring in the shrimp ponds more than the irrigating tidal channels, the amount of evaporation has little overall influence on major element concentrations in solution (Additional file 1: Fig. A5). Evaporation appears low in May (≤ 10%) and does not significantly affect element concentrations in shrimp ponds. This is likely a result of regular shrimp pond drainage of wastes and replenishment from new tidal channel water, which would reduce water residence time and the effects of evaporation. However, the evaporation models were based on a sample close to high tide and are thus minimum estimates of evaporation, although in the nearby Bangladesh Sundarbans, salinity changes between high and low tide in the dry season tidal channels are often around ~ 1ppt or less, with some tidal channels seeing practically no salinity change between tides at all [53,54].
Trace elements of interest, such as As and Se, also do not show any clear trends with stable isotopes (Additional file 1: Fig. A10). Weak to moderate relationships are seen with all other elements as well (Additional file 1: Fig. A11). Essentially, evaporation is not great enough to overcome the large variance in starting compositions of the shrimp pond water, or the effects of nonconservative element behavior are larger than the effects of evaporation, effectively erasing any possible correlation between stable isotopes and nonconservative behaving elements.
DOC shows a strong positive relationship with stable isotopes (Fig. 7). As δ 18 O and δ 2 H isotopes become heavier in shrimp ponds, DOC increases as well, likely because these variables are affected by endogenous processes in shrimp ponds (e.g., biological processes or evaporation), while other geochemical parameters are more affected by exogenous processes (e.g., upstream rock weathering, seawater mixing). For example, as shrimp pond waters undergo more evaporation, they become more stagnant with lower flow conditions, algaerich, and produce DOC through algal photosynthetic activity (e.g., [19]). Additionally, all shrimp ponds have similar DOC concentrations in their irrigation source (as seen with little variation in tidal channel DOC values in Fig. 7), which allows identification of the small changes caused by evaporation. DOC values do not show a relationship with isotopes in tidal channels because tidal channel isotopic values are predominantly influenced by rainwater-seawater mixing instead of evaporation.
Lastly, it must be noted that expanding isotopic analysis of shrimp pond waters to other months/seasons may be useful for helping determine geochemical relationships as well, particularly because when excluding isotopes, multidimensional scaling (MDS) shows July shrimp ponds are distinct from other samples (Additional file 1: Fig. A12).

Multiple linear regression and predictors of As and Se
Nearly all our samples plotted similarly to seawater in a Piper diagram (Additional file 1: Fig. A13), and were thus collectively used in multiple linear regression to assess important predictor variables for As and Se, even though most elements showed large variability in concentration.
Multiple linear regression illustrates that Ni, P and V are important predictor variables for Se and/or As in most water samples (Fig. 8). Phosphorus is regarded as a "chemical analog" to As, and like As (as pentavalent arsenate) in oxidized aqueous environments, P is mainly in an oxyanion form (as PO 4 3− ) [69]. However, in the main predictive model for As, P is positively correlated (Additional file 1: Table A3), which would suggest that even though both elements are likely negatively charged species in solution, their concentrations are not explained by competitive adsorption. Vanadium is likely predominantly an oxyanion as well given the pH and oxidizing conditions of the surface water samples, with the most common ion of V being H 2 VO 4 − in natural waters [20]. Because V is positively correlated with As in the predictive model (Additional file 1: Table A3), it may also behave similarly in solution but not compete directly for sorption sites. However, Ni is most often a cation species with a valence of + 2 in natural waters [39]. When analyzing primary influencing variables for As (Cu, V, Ni and P), it is clear that only Ni is negatively correlated, illustrating that Ni may improve the model because of different chemical behaviors (i.e., sorption affinities) as a cation in solution (Additional file 1: Table A3).
Selenium also exists as oxyanions in oxidizing environments such as surface waters [73]. When looking at the predictive variables V, Ni, and DOC for Se,V has a negative coefficient in the model while Ni has a positive coefficient (Additional file 1: Table A5). This is opposite that of the main As model (Additional file 1: Table A3). Thus, for Se, even though V is likely an anion in solution and Ni a cation, the anion is associated with a decrease in Se in the model. This may be due to competitive sorption properties of Se and V on solid particle surfaces (particularly with selenite on Fe-oxyhydroxides or clay minerals), which would increase V in solution as Se sorbs. Ni may complex with aqueous species that behave similarly to Se but do not compete with sorption sites.
pH is important to include in modeling because of its known effects on speciation and mobility of Se and As (e.g., [69,73]). However, models that include pH do not significantly improve As or Se prediction (Fig. 8), potentially due to relatively small variance in pH values (Table 1). Thus, pH does not likely have a large, statistically significant influence on either Se or As concentrations in these surface waters, although it is acknowledged that pH is important in speciation of these elements and thus their chemical characteristics in solution.
When adding isotopes to predictive modeling to determine whether evaporation/mixing may affect Se and As predictability, the adjusted R 2 values for both Se and As increase to both 0.70 (p = 0.016) and 0.90 (p = 0.0041), respectively (Additional file 1: Tables A6, A7), but only 14 samples are included in the analysis, omitting all July shrimp ponds. The influence of isotopes must be taken with care, as clear bivariate relationships between Se, As and δ 18 O (and δ 2 H) are not apparent (Additional file 1: Fig. A10), and the inclusion of both isotopes (δ 18 O and δ 2 H) in the Se model (Additional file 1: Table A7) clearly violate the multiple linear regression assumption of no multicollinearity among predictor variables. However, future work investigating the predictive nature of stable isotopes for trace elements in surface water or groundwater modeling may be warranted because as proxies of mixing or evaporation, O and H isotopes may provide additional insight as to where certain elements originated from or if evaporation concentrated elements in solution.
Additionally, because organic matter such as DOC can affect trace elements in multiple ways, such as complexing with As species or competing for sorption sites (e.g., [41,72]), it was another geochemical parameter worth examining. For Se, it was found that the predictor variables DOC, V, and Ni resulted in an adjusted R 2 of 0.67 (p < 0.0001) (Additional file 1: Table A5). Thus, Se may also be affected by DOC, such as by surface water reoxidation of biogenic elemental Se (BioSe) that originated from microbes [73] in high DOC waters. More research in shrimp pond elemental cycling would help shed insight on this. Lastly, it is apparent that adding Cl as a representative saltwater ion to the variables DOC, V and Ni improves the predictive model (Additional file 1: Table A4), which may be indicative of slight seawater mixing or seasonality effects on aqueous Se. However, it is noted that Cl does have a strong positive correlation with V (Additional file 1: Fig. A1), and thus violates the assumption of no multicollinearity for multiple linear regression.

Applicability of Se and As models
While these models do not identify the mechanisms that control As and Se concentrations in shrimp ponds and tidal channel waters, they are very useful in gauging what elements/geochemical parameters are most important in influencing As and Se concentrations in surface waters. This extends to upstream (non-headwaters) Ganges River samples [15], their Table S4), where Cu, Vi, and Ni result in a great multiple linear predictive model fit for As (adjusted R 2 = 0.64, p = 7.2e-15) (Additional file 1: Table A8). Additionally, when examining tidal channel samples from Ayers et al. [8], predicting As with Cu, V, Ni, and P results in a great overall fit as well (adjusted R 2 = 0 0.74, p = 0.024) (Additional file 1: Table A9). However, Se predictive modeling for tidal channels in Ayers et al. [8] is much poorer, with DOC, V, and Ni only resulting in an adjusted R 2 of 0.24 (p = 0.17). Nevertheless, the general approach of multiple linear regression with several important predictor variables such as Ni and V or other geochemical analogues appears to be applicable to other surface water bodies, particularly rivers, and can be especially useful in better understanding and predicting the concentration of As in solution. Further research in quantifying Se in surface waters will help in understanding what controls/predicts its geochemical variability.
These models also provide a way to estimate As and Se concentrations in Bangladesh surface waters without direct measurements of As or Se. Model validation using other published datasets indicates that the model parameters obtained in this study do not produce model predictions that agree well with observed As and Se concentrations (i.e., Additional file 1: Fig. A14). However, the general approach used here can help improve predictions and measurements in other regions throughout the world, particularly when several geochemical parameters may be missing from samples within an area. For example, if a study has only analyzed As or Se in 15% of samples, but the remaining samples contain potentially important predictor variables such as Ni or V, a modeling solution can be applied to estimate As or Se in the remaining samples.
Although this study assumed linear relationships in predictive modeling, future work involving nonlinear predictive modeling techniques such as machine learning techniques (i.e., boosted regression trees) is warranted, as these techniques have been shown to outperform multiple linear regression predictive modeling in groundwater [49].

Selenium and arsenic antagonistic relationship
Se is known to combat the toxic effects of As, and As is known to inhibit Se toxicity (e.g., [30,67,68]). Elevated Se and As can be taken up by biological organisms, and may eventually sorb onto soils/sediment. This is particularly relevant for As in the area, as many shrimp ponds are utilized for rice farming in the wet season [6], where flooded soils could remobilize sorbed As through reductive dissolution of Fe oxyhydroxides (e.g., [57]). However, even if Se becomes elevated in soils from tidal channel irrigation, the reducing conditions in flooded rice paddy soils could lead to immobile elemental selenium being formed, which is not bioavailable [73]. Thus, rice could potentially take up arsenic, but not selenium, with humans therefore lacking the detoxifying effects of selenium when ingesting the rice. Further work is needed to explore the cycling and antagonistic relationship between Se and As in Southwest Bangladesh, such as identifying if sediments in contact with Se and As-rich waters are elevated in both those elements, what fraction is bioavailable, and if Se concentrations in tidal channel waters in the peak monsoon season become heavily diluted. Although much research in Bangladesh has examined As in crops such as rice (e.g., [37,71,72]), crop/aquaculture Se should be more heavily researched in Southwest Bangladesh as well, particularly where Se water concentrations are high.

Conclusion
Through examining surface water chemistry in Southwest Bangladesh in two different months over a wide spatial area, it was deduced that: (1) Monthly precipitation differences between May and July has a greater effect on tidal channel water composition compared to shrimp pond composition; (2) There is a large compositional difference between shrimp ponds and the tidal channels irrigating them in the early monsoon (July), but not in May; (3) Evaporation in shrimp ponds and endogenous pond effects in general have a relatively minimal impact on surface water trace element and major ion chemistry based on May pond data, although DOC has a strong positive correlation with δ 18 O and δ 2 H isotopes and thus evaporation; (4) Arsenic and selenium concentrations are elevated above WHO drinking water guidelines in the majority of surface water samples and are correlated with other trace elements in solution such as Ni and V; and (5) Predictive modeling of hazardous trace elements in surface waters may prove useful in future studies throughout the world when measurements of certain toxic elements cannot always be easily made.