Physical controls and ENSO event influence on weathering in the Panama Canal Watershed

Recent empirical studies have documented the importance of tropical mountainous rivers on global silicate weathering and suspended sediment transport. Such field studies are typically based on limited temporal data, leaving uncertainty in the strength of observed relationships with controlling parameters over the long term. A deficiency of long-term data also prevents determination of the impact that multi-year or decadal climate patterns, such as the El Niño Southern Oscillation (ENSO), might have on weathering fluxes. Here we analyze an 18-year hydrochemical dataset for eight sub-basins of the Panama Canal Watershed of high-temporal frequency collected between 1998 and 2015 to address these knowledge gaps. We identified a strongly positive covariance of both cation (Ca2+, Mg2+, K+, Na+) and suspended sediment yields with precipitation and extent of forest cover, whereas we observed negative relationships with temperature and mosaic landcover. We also confirmed a statistical relationship between seasonality, ENSO, and river discharge, with significantly higher values occurring during La Niña events. These findings emphasize the importance that long-term datasets have on identifying short-term influences on chemical and physical weathering rates, especially, in ENSO-influenced regions.

Empirical studies over the past two decades have documented the importance of small mountainous rivers (SMRs), particularly those in the tropics, on global silicate weathering and associated CO 2 consumption budgets. In a search for controls, studies have recognized the strong positive feedback between physical and chemical weathering [1][2][3] and noted the importance of underlying volcanic lithologies on maintaining elevated chemical yields [4][5][6][7][8][9][10][11] . Others have identified strong correlations/associations between weathering yields and precipitation/ runoff 2,3,9,10,12 and, more recently, land use/landcover practices 10,11,13 . Yet, despite their important contribution to silicate weathering, few datasets for SMRs exist at high temporal resolution 10,14 and/or duration 15 for such high-yield terrains. This lack of long-term records of high temporal resolution requires caution in generalizing the aforementioned statistical relationships, as years of anomalously high or low precipitation can result in misinterpretation of silicate weathering rates. This is particularly important as SMRs are presently being incorporated into newer global chemical weathering and associated CO 2 drawdown models [16][17][18] . Furthermore, the absence of reliable long-term datasets prevents insight into the impact that multi-year or decadal climate patterns, such as the Pacific Decadal Oscillation (PDO) and El Niño Southern Oscillation (ENSO), can have on weathering fluxes.
In the tropics of the Caribbean and South America, ENSO events have been shown to have a pronounced and varied effect on rainfall [19][20][21][22] , streamflow [22][23][24] , soil moisture 23 , and vegetation 23,25 . For example, a review of hydrological data for Colombia observed that ENSO effects between 1997 and 1999 were stronger for stream flow than for precipitation, due to concomitant effects on soil moisture and evapotranspiration, with lower than normal soil moisture and stream flow occurring during El Niño conditions and the reverse situation pertaining during La Niña episodes 23 . Similar trends with ENSO events have been observed in Panama. For example, a previous study documented that an average of 8% less rainfall was received across almost all regions of Panama during 13 El Niño episodes between from 1920 to 1983 19 . Such ENSO-related reductions in precipitation and consequent runoff [26][27][28] 29 . It follows that these pronounced hydrological anomalies should also impact weathering fluxes, but this has not been documented to date. The Panama Canal Watershed (PCW; 2,982 km 2 ; 9° N, 80° W), comprised predominantly of eight major subwatersheds, offers an ideal location to evaluate the physical and climatic controls on long-term weathering rates and CO 2 consumption in tropical SMRs (Fig. 1). For example, the largely forested, steeply-sloping sub-watersheds on the north side of the canal-Río Gatún, Río Boquerón, Río Pequení, Río Chagres, and Río Indio Este, differ markedly from their mostly deforested, gently-sloping counterparts to the south-Río Cano Quebrado, Río Trinidad, and Río Cirí Grande 30,31 . The region also exhibits a strong trans-isthmus rainfall gradient, with the amount of precipitation received on the windward Atlantic coast more than twice that received on the leeward Pacific coast (~ 4,000 mm/year versus < 1,800 mm/year) 32 . Lastly, pronounced differences in regional precipitation values during ENSO events permit an evaluation of the impact of short-term climate patterns on weathering fluxes 21,26 .
Here we utilize a robust, high-temporal frequency, hydrochemical dataset collected over 18 years between 1998 and 2015 by the Panama Canal Authority (ACP) to calculate annual and long-term cation and sediment fluxes for the PCW. We compare the resulting fluxes to potential controlling variables, such as mean annual rainfall and  www.nature.com/scientificreports/ ering rates. For example, a two-tailed t-test revealed significant differences in both cation (α = 0.05, p = 0.038) and suspended sediment (α = 0.05, p = 0.011) yields between the north and south watersheds despite no noted differences in average discharge values (α = 0.05, p = 0.28). Our observed lack of statistically significant differences in discharge points to the importance of lithology, as well as other variables, in maintaining high chemical erosion rates. For example, watersheds on the north side of the canal are largely underlain by mafic to intermediate volcanic rocks compared to the largely sedimentary cover for those on the south side 11 . The importance of volcanic lithology in maintaining disproportionately high cation yields was previously documented as part of an isthmus-wide study 10,11 . While we observed a negative, and somewhat variable statistically significant relationship between cation yields and basin-wide mean annual temperature (r avg = − 0.39), this is likely due to increased cloud cover associated with precipitation events. This idea is supported by an observed negative statistical relationship between basin wide mean annual rainfall and temperature (r = − 0.47, p < 0.05) over the study period. Previous studies have documented a strong relationship between physical and chemical weathering rates worldwide in both large catchments 33 and SMR watersheds [1][2][3]10 . Analyzed collectively, our watersheds confirmed this linkage, with 7 of the 11 study years exhibiting a significant positive relationship between cation and sediment fluxes (r avg = 0.62, p avg ≤ 0.17) (Supplementary Tables 5 and 6). However, we observed notable differences in correlation strength when watersheds were analyzed individually over time. For example, we observed positive, yet slightly weaker correlations between these two parameters in the Chagres (r = 0.69, p = 0.02) and Boqueron (r = 0.69, p = 0.02) watersheds despite their location on the windward side of Atlantic range where the precipitation rate is corresponding high (> 3,000 mm/year) (Supplementary Table 4) This counterintuitive result could be explained by the fact that these watersheds are characterized by heavy forest cover (> 90%), which has been previously shown to increase infiltration and soil strength and, therefore, reduce surface runoff in tropical watersheds 34,35 .

Relationship of weathering rates with land use/landcover (LULC) practices. Comparisons
between weathering fluxes and LULC practices in tropical SMR watersheds have been limited in scale to date 10,11,13 , due to low temporal resolution geospatial data. However, annual LULC data [forest and mosaic (i.e., forest plus croplands)] measured for the PCW on an annual scale by the ESA Climate Change Initiative (CCI) allowed for their direct comparison with cation values over the 17-year period (The presence of large-scale cropland in only one of the watersheds prevented a comparison with LULC practice). Forest cover was the dominant LULC (> 50%) in all but one of our watersheds and comprised upwards of 95% of LULC in three watersheds (Supplementary Table 7). Our Pearson statistical analysis revealed a positive relationship (r ≥ 0.47, r avg = 0.63) between percent forest cover and cation fluxes for the entire study period (Supplementary Table 5). Alternatively, we observed mosaic LULC practices to exhibit a strong statistically-significant negative relationship (r ≤ − 0.54, r avg = − 0.74) with cation fluxes, with 11 of 18 years exhibiting statistical significance.
Our findings are in general agreement with those from a previous isthmus-wide study for Panama based on spot sampling or rivers over a decade 10,11 and suggests land use practices have fundamentally altered hydrological flow pathways and water residence times in this tropical hydrological system. Previous studies in tropical forested watersheds have observed that some ~ 80-90% of rainfall infiltrates into the soil 34 , where roots can retain water during periods of soil saturation and release it throughout the dry season as baseflow, thus reducing total forest catchment runoff, but increasing flow consistency 34,35 . The concomitant increase in water residence time in contact with fresh mineral surfaces in the soil regolith allows for increased mineral dissolution and solute export. Conversely, soil compaction in agricultural pasture and croplands increase the likelihood of runoff, which is supported by previous studies in the PCW that show mosaic catchments produce 1.8 times more runoff than their forested counterparts 35,36 . Interestingly, the decrease in correlation strength between amount of forest cover and cation fluxes in our study watersheds over time coincided with an increase in percent forest cover in the study catchments over the same period ( Fig. 3; Supplementary Tables 5 and 8). This counterintuitive pattern may be attributed to the fact that trees in abandoned mosaic plots, reverting back to forest, have not yet established root networks capable of altering hydrological flow pathways. Impact of ENSO on stream discharge and weathering fluxes. We hypothesized that both seasonality and ENSO conditions would influence stream discharge (and by cation and sediment fluxes) and, furthermore, that the two factors would interact. To evaluate the effect of ENSO conditions on stream discharge and weathering fluxes, we used the Oceanic Niño Index (ONI), defined as a 3-month running mean of SST anomalies in the Niño 3.4 region (5º S-58º N, 170º W-120º W) 37 . We justify the use of ONI versus the Southern Oscillation Index (SOI), as previous global compilations studies have observed good agreement between these two indices 28 and that ONI allows for distinct classification of ENSO events (The National Oceanic and Atmospheric Association classifies El Niño and La Niña events as 3-month ONI running means either < − 0.5 or > + 0.5 than the long-term average, respectively). We modelled the influence of seasonal periods and ONI on river-discharge and weathering fluxes for the pooled dataset using a mixed-effects model, with seasonal interval (i.e., JFM, etc.) and ONI as fixed effects, and river (i.e., watershed) as a random effect. We accounted for temporal auto-correlation inherent in time-series sampling using an autoregressive AR(1) correlation-structure for residuals. Our model identified a significant interaction of seasonality and ONI on river discharge (p < 0.0001) and cation fluxes (p < 0.0001) and near significance on sediment fluxes (p = 0.056) (Supplementary Table 9). This decoupling between river discharge and sediment flux implies that discharge is not the sole process controlling the sediment flux of PCW rivers. A positive slope-estimate for seasonality is consistent with/reflects the general increase in precipitation from the dry to wet season transition, and corresponding increase in discharge. A negative slope-estimate for ONI suggests less discharge with higher ONI values, which agrees with previous hydrological analyses for the region 23 www.nature.com/scientificreports/ Our visual inspection of the PCW dataset identified elevated mean discharge and weathering fluxes for several seasonal timesteps during La Nina conditions (Fig. 4). We then performed a series of one-way ANOVA tests to determine statistically significant differences (p < 0.05) between La Niña, El Niño and neutral discharge as well as weathering fluxes. We confirmed statistical differences only for the December-January-February (DJF) and January-February-March (JFM) periods; with La Niña events exhibiting the highest average values for 9 of the 12 tri-monthly time periods (Fig. 4). Our subsequent series of ad-hoc Tukey tests on the dataset revealed significant differences between La Niña and El Niño (p < 0.001) and La Niña and neutral (p ≤ 0.001) discharge as well as weathering flux values for the DJF period (Fig. 5). Significant differences were also observed between weathering flux values and La Niña and neutral discharge (p ≤ 0.035) and La Nina and El Niño (p ≤ 0.002) discharge values for the JFM period. Previous hydrological studies in Panama and northern South America have identified a positive correlation between monthly average and/or daily maximum precipitation 23,38 , and stream discharge values 23,28,38,39 with DJF SOI (i.e. drier El Niño and wetter La Niña); the time period when El Niño events typically achieve their maximum. In neighboring Columbia, strong seasonal moisture advection anomalies created by the winds of the "CHOCO jet" have been offered as a possible explanation for a robust regional relationship between precipitation and the DJF SOI time interval 23 , while weaker correlations for March-April-May are potentially explained by the fact that ENSO is either just starting to develop or is declining at that time of the year 23 .
Our observed ENSO-driven deviations in stream discharge have a demonstrable impact on annual cation weathering fluxes. For example, the four highest basin-wide annual average cation yields (Supplementary Table 1; watersheds with 18 year record, only) were recorded in years dominated by strong La Niña events (1999, 2007, 2010, and 2011). This trend was also supported by the sediment yield records for overlapping years (2007, 2010, 2012; Table S2). We also observed that the timing and relative strength of the La Niña event is also an important factor as these same four years exhibited ONI values < 1 through most of the wet season. The only exception to this pattern was in 1998, which was also marked by a strong El Niño event (ONI values > 1) during the first half of the year. Furthermore, a particularly strong La Niña event that occurred during the entirety of the 2010 wet season resulted in average annual cation and sediment fluxes (15-76% and 18%-404%, respectively) that were substantially greater than their respective long-term averages. Exceptionally high rainfall during a 7-8 December 2010 storm, having a return period estimate of 2,000 year 40 , resulted in numerous landslides throughout the PCW 40,41 and river sediment loads that overwhelmed the Panama City water treatment facilities 40 . High rainfall events associated with La Niña periods have also been linked to an increase in landslides in Columbia and Venezuela 42,43 . Landslides have been previously shown to play a critical role in maintaining the disproportionately high sediment and cation yields in SMR catchments 44,45 , and thus are likely playing a similar role in the PCW.
Conversely, the strong El Niño event observed throughout the 2015 calendar year coincided with our lowest basin-wide average annual cation weathering rate over the 18-year study period. This event was also marked by anomalously high ONI values (i.e. > 2) during the latter half of the year and anecdotal stories regarding water shortages across the canal zone. However, we observed a more variable hydrological response to both neutral and El Niño conditions/events throughout the remainder of the study period, which is supported by our statistical analysis. Previous studies have documented an 8% decrease in rainfall in almost all regions of Panama during www.nature.com/scientificreports/ 13 El Niño episodes between 1920 and 1983 19 and below median Lake Gatun inflow in 17 of 20 instances when SST anomalies in the NINO3 region exceed 0.6 ºC 27 . While we did not observe a similar overall hydrological response to El Niño events, this might be in part due to the relative lack of strong El Niño events over our 18-year study period 37 . ENSO driven deviations in weathering rates are also apparent through comparisons to previous studies utilizing shorter-term data sets. For example, a recent determination of cation weathering rates for the Chagres and Pequini watersheds 10,11 calculated long-term values 80% and 48% greater, than those of this study. Unsurprisingly, 63% of the spot samples and instantaneous discharge values we used to construct the associated weathering equations for the study were collected during La Niña events. With atmospheric modeling predicting an increase in the frequency of both extreme El Niño 46 and La Niña events 47 due to greenhouse gas warming, caution will need to be employed when evaluating future empirically-based weathering studies.
The new findings presented here not only confirm the need for long-term weathering studies in tropical regions, but also suggest that caution should be employed when incorporating data from regions influenced by multi-year or decadal climate patterns into global compilation studies. While much progress has been made over recent decades on the determination of weathering fluxes from high-yielding terrains in the Caribbean 7,9,15 , Central America 10,11 , southeast Asia 8 , and Oceania 2 , these are all regions heavily influenced by ENSO events. Interestingly, these regions have also been shown to play a disproportionate role in the delivery of dissolved 48,49 and particulate organic carbon 48 to the global ocean, and by extension, other nutrients such as PO 4 . While ENSO events have been shown to affect vertical mixing and associated upwelling of nutrients in affected ocean areas 50 , our data suggests nearshore locales will also be impacted by concomitant changes in nutrient delivery from the coast. For example, decreased upwelling of nutrients in the eastern Pacific during El Niño periods would be compounded with a corresponding decrease in nutrient export from land, thus exacerbating nutrient limitation in these locales.
Prior to its use in denudation calculations, streamwater cation concentrations were corrected for sea-salt contribution as follows: non-sea-salt concentration = measured concentration -(sea-salt correction ratio to www.nature.com/scientificreports/ Cl -)*(Cl -). Following Murphy and Stallard (2012), we used the following species-to-chloride ratios for this adjustment: Na + /Cl -= 0.85251, K + /Cl -= 0.01790, Mg 2+ /Cl -= 0.09689, and Ca 2+ /Cl -= 0.01879. We further adjusted for non-silicate contributions of Ca and Mg using volcanic end-member ratios of 0.5 for Ca/Na and 0.5 for Mg/ Na, previously established 33 . Our use of these values instead of a continental granitic counterpart is justified by the high Ca/Na and Mg/Na concentrations in both waters and the mafic to andesitic character of igneous rocks across central Panama 10,11 . Using previously established methodology 10 , we employed a multistep process whereby individual cation concentrations in the data set were first multiplied by the corresponding average daily discharge value in order to produce an instantaneous chemical denudation value. We then prepared x-y plots of instantaneous denudation values with respective discharge to produce specific elemental yield determination equations (Supplementary  Table 10). Our approach is supported by high average correlation (r 2 ) values observed for the instantaneous denudation value-discharge comparisons (Ca = 0.80; Mg = 0.85; Na = 0.82; and K = 0.89). Finally, we substituted daily discharge values over each of the 18 years of record into the equations, and our calculated denudation values were subsequently divided by watershed area to produce annual and long-term estimates of cation weathering rates. For the suspended sediment calculations, we multiplied average daily discharge values by the corresponding average daily suspended sediment concentrations to produce a daily sediment load. Our daily suspended sediment loads were then compiled and subsequently divided by watershed area to produce annual and long-term estimates of suspended sediment yields. The ALOS Global Digital Surface Model for the Greater Panama Canal Watershed study area was hydrologically conditioned to reduce noise 56 . Multiple flow direction (MFD) flow accumulation was computed over the hydrologically conditioned digital surface model (DSM) using an A T least-cost search algorithm to traverse depressions and obtacles 57 . We then extracted the stream network from the DSM and flow accumulation, and the stream gage stations were snapped onto the stream network. We derived watershed basin outlets at the stream gages from the flow direction of the stream network. Landcover data from the ESA CCI Landcover Dataset for 1998 through 2015 were reclassified and re-categorized as shown in Supplementary Tables 11 and 12. We computed topographic, hydrological, and landcover analyses for each basin. The topographic parameters analyzed included elevation, slope, and aspect. Our hydrological parameters included flow accumulation, stream geometry, stream distance, stream order, and stream statistics. Our landcover parameter was the percentage of landcover in each class over the study period. Our topographic, hydrological, and landcover maps were visualized with shaded relief based on the composite of direct illumination derived from topographic relief and diffuse illumination derived from the sky-view factor 58 .

Statistical analysis.
We performed all statistical tests using JMP (Pro version 14.0.0; SAS Institute, Cary, NC). Data which did not meet conditions for normality using a Shapiro-Wilk check test were log transformed prior to analysis. Sample size (n) associated with our statistical analyses are provided in the respective supplementary data tables, unless otherwise noted.
We modelled the influence of seasonal period (average tri-monthly interval) and ONI 37 on the river-discharge (pooled tri-monthly measures n = 1,284), cation flux (n = 1,284), and sediment flux (n = 780) datasets using a linear mixed-effects model (Eqs. 1, 2 and 3): where Q tr, , Cat tr , and Sed tr , represents the estimated discharge, cation flux, and sediment flux, respectively, for a particular tri-monthly period over the 18-year dataset; t, for a particular river; r (Only streams with a full 18 years of data were included in our statistical model for discharge and cation fluxes, which resulted in a total of six river basis. The included sediment record spans eleven years for the same basins). These discharge and weathering flux values were modelled as a function of TMI, a nominal variable indicating the tri-monthly seasonal interval within a year, ONI, a continuous variable, and ONI*TMI representing their interaction. Both visual inspection (graph boxplots of ONI as y-axis and DJF, etc. as x-axis, for the 18-year dataset) and Kendall's tau suggested little association between ONI and trimonthly interval (T b = − 0.0385, p = 0.43), therefore we used both as explanatory variables in our statistical analyses. In addition to these fixed effects, µr represents the random effect due to discharge values being associated with the rth river and ε is the random error. River was designated as a random factor, as they are a small sample of the total number found in Panama and were not chosen with some explicit comparison in mind.