2.1. Site description
La Selva Biological Station, Costa Rica (LSBS), is a 1600 ha tropical wet forest reserve, with elevation ranging from 22 to 146 meters above sea level (Fig 1). LSBS is located at the transition from the upland foothills of the Cordillera Central of Costa Rica and the Caribbean lowlands, and is the lowland terminus of the altitudinal transect in the Braulio Carrillo National Park (Fig 1). Mean annual rainfall from is 4300 mm, with dry season January to April and wet season July to December (Sanford et al., 1994).
We focused on the lower 75 m reach of the Taconazo stream (10.432, -84.013) to partition and quantify CO2 input (upstream input, NEP, and GWCO2) and output (FCO2, downstream hydrologic export) fluxes. The Taconazo watershed area is 0.28 km2, all of which is forested with a closed canopy (Table 1). The Taconazo is a low solute stream at LSBS, with long-term (1997 - 2015) mean water temperature of 25.1 ˚C, pH of 5.5, discharge of 0.06 m3 s-1, and NO3- and SRP of 192.4 µg L-1 and 4.9 µg L-1, respectively (Ganong et al. 2015). The Taconazo has received significant study as a model headwater lowland tropical stream. Numerous studies from the Taconazo have contributed to better understand of the hydrology (Genereux et al. 2005; Genereux and Jordan 2006), biogeochemistry (Small et al. 2012; Oviedo-Vargas et al. 2015; Osburn et al. 2018), and ecology (Ardón et al., 2006; Ramírez et al., 2003; Rosemond et al., 2002) of tropical lowland headwater streams.
2.2. Data collection
We collected stream data from April 1, 2013 - September 30, 2013. Discharge (m3 s-1) was continuously measured from the V-notch weir constructed near the outflow of the Taconazo; further descriptions of the weir in Zanon et al. (2014). Rainfall (mm) was collected from the LSBS weather station, located ~900 m from the Taconazo. Meteorological data are available at: https://archive.tropicalstudies.org/meteoro/default.php?pestacion=201.
Groundwater discharge into the 75 m study reach was measured using conservative tracer injections to the stream. In the dry season (April), groundwater flux into the reach was 20% of total discharge at the end of the reach (Oviedo-Vargas et al. 2015), compared to the wet season (May-September) flux of 9% (Oviedo-Vargas et al. 2015, Ardón et al. 2013). For the dry and wet season, we multiplied discharge by the respective percentages as an estimate of groundwater discharge (QGW). To calculate discharge into the upstream end of the study reach (Qu), the groundwater discharge in the reach was subtracted from stream discharge measured at the end of the study reach (Qd).
A YSI 600xlm multiprobe sonde was secured 5 cm above the stream bottom in a PVC tube with holes drilled to allow for exchange of stream water. The multiprobe continuously measured temperature (˚C), pH, and DO (mg L-1 and %-saturation) at hourly intervals. The sonde was retrieved every two weeks for recalibration and to clean fouling on the sensor heads. Partial pressure of CO2 (ppmv) was measured using a Vaisala GMT221 infrared gas analyzer (IRGA) in the stream and in a riparian well. The sensors were prepared for submerged deployment as described in Johnson et al. (2010) and pCO2 was logged hourly using a Campbell CR1000 datalogger. pCO2 data were processed according to Johnson et al. (2010) for depth and temperature. We removed data from all sensors when the water level fell below the height of the sonde (e.g. conductivity < 0.005, DO near air saturation). Sensor data underwent further QAQC following the guidance of Taylor and Loescher (2013) and using the datacleanr R package (Hurley 2021).
2.3. Evaluating aqueous concentrations of CO2 and O2
We evaluated the departure of stream CO2 and O2 concentrations relative to concentrations expected at atmospheric equilibrium, and the magnitude of gas departure over time. Departure of O2 and CO2 from equilibrium with the atmosphere in freshwaters is mediated by biochemical reactions that directly link O2 and CO2 (e.g. aerobic metabolism), water temperature, bicarbonate chemistry, and gas exchange (Vachon et al. 2020). We calculated CO2 departure as the difference of measured aqueous concentration ([CO2]aq) minus the concentration at atmospheric equilibrium at 400 ppmv ([CO2]sat) using a temperature-corrected Henry’s Law constant. O2 departure was calculated as the difference between in-stream concentration minus the concentration at saturation which is calculated as a function of temperature and barometric pressure, following the outline in Hall and Hotchkiss (2017) and detailed below.
Paired O2 and CO2 departure were plotted as a data cloud, O2 departure vs. CO2 departure. To evaluate movement of the cloud over time, we estimated 95% confidence interval ellipses for each month of data collection using the ellipse v0.4.2 R package (Murdoch and Chow 2020). In each ellipse, we calculated the: 1) centroid, which indicates the location in Cartesian space for both O2 and CO2; 2) 1/[slope] through each ellipse, which informs the efficiency of metabolism between O2 and CO2, interpreted as the moles of CO2 produced per moles of O2 consumed through ecosystem respiration. See Vachon et al. (2020) for more details on these metrics.
To examine the role of individual rainfall events (2067 mm over 180 days) on pCO2, we selected three storm events: one event from the dry season (April), early wet season (June), and late wet season (August). For the three events, we examined the hysteresis response of pCO2 to discharge. For each event, we selected data 1 hour prior to the rising limb, through the peak in discharge, and for 5 hours following return to pre-storm discharge and before the next rising limb in the hydrograph (Evans and Davies 1998; Dinsmore and Billett 2008). We evaluated each hysteresis response for diagnostics describing the solute-discharge rotation (clockwise vs. counterclockwise vs. ‘figure 8’) and trend (positive vs. negative) (Evans and Davies 1998).
2.4. CO2 flux calculations and reach scale DIC mass balance
Our data collection allowed for independent calculation of CO2 input and output fluxes in the study reach. We calculated aqueous CO2 from the in-stream (see section above) and riparian well ([CO2]GW) stations from pCO2 data using the temperature dependent Henry’s Law constant from Plummer and Busenberg (1982). Stream CO2 at saturation ([CO2]sat) was determined using the same method, taking 400 ppm as the atmospheric CO2 concentration (Rocher‐Ros et al. 2019).
We compared sensor estimates of DIC to DIC estimates from samples of stream water. Samples were collected in 10 mL syringes and 4 mL of water was injected into an inverted sealed serum vial with a 22-gauge needle and 0.45 µm pore filter. Samples were equilibrated on a shaker table, and 250 µL headspace was removed and pCO2 analyzed on an SRI Instruments gas chromatograph (Las Vegas, Nevada, USA). Conversions from pCO2 to DIC were calculated as described in the supplementary file for sensor data. We matched the dates of syringe samples with sensor derived estimates and compared the two methods using the Student’s t-test.
We defined CO2 inputs to the study reach as the groundwater CO2 flux (GWCO2) and in-stream net ecosystem production (NEP). Groundwater CO2 flux was calculated as
Where SA is the benthic area of the stream and the flux units are mol C m-2 h-1
CO2 inputs from in-stream aerobic metabolism was estimated as hourly NEP using the instantaneous direct metabolism method outlined in Hall and Hotchkiss (2017). Briefly, NEP was estimated from the volumetric O2 mass balance equation:
Where NEPi is hourly instantaneous metabolism (mol O2 m-2 h-1), Oi is the O2 concentration at each timepoint (mol m-3), Osat,i is the concentration of O2 at saturation for the same timepoint (mol m-3), z is stream depth (m), and KO is the gas exchange coefficient for O2 (h-1). We converted KO from KCO2 (see equation 3) via Schmidt scaling (Raymond et al. 2012). NEP was converted to mol C m-2 h-1 assuming a 1:1 respiratory quotient between O2 and CO2 (Rocher-Ros et al. 2019). We selected the direct metabolism method for determining NEP because the lack of a diel O2 signal (Fig S1) and relatively high reaeration make the Taconazo ill-suited to alternative stream metabolism methods (Appling et al. 2018). Further, the direct method allows for estimation of NEP at finer temporal compared to alternative approaches, which are resolved at the daily scale.
CO2 output as areal CO2 evasion (FCO2, mol C per stream surface area per hour, m-2 h-1) was calculated as
Where KCO2 is the dry (1.13 h-1) or wet season (0.43 h-1) gas exchange coefficient reported in Oviedo-Vargas et al. (2015), and is depth in meters measured at the weir.
We conducted a volumetric mass balance in terms of DIC at the reach scale to compare reach scale estimates with watershed scale hydrologic inputs (upstream inputs and downstream export). The mass balance scales the fluxes calculated in equations 1-3 from mol C m-2 h-1 to reach as mol DIC h-1; conversion of CO2 to DIC are detailed in the supplementary file. Upstream DIC inputs to the study reach were calculated as the product of upstream discharge and [DIC] and upstream discharge was as the difference of downstream discharge and groundwater flux in the reach, which is an assumption supported in previous work (Ardón et al. 2013, Oviedo-Vargas et al. 2015). Thus, the full mass balance is:
Where C is DIC (mol m-3), Q is discharge (m3 h-1), L is the reach length (75 m), w is wetted width (1.5 m in the dry season, 2.8 m in the wet season), z is depth (m) measured at the weir, NEP is volumetric net ecosystem production (mol DIC m-3 h-1) and subscripts d, u, and GW refer to downstream, upstream, and groundwater, respectively. We evaluated the mass balance through the dDIC dt-1 term. If all relevant fluxes are accounted for, the change in should average zero, representing DIC inputs and outputs in balance with no net change in C storage in the reach.
We aggregated fluxes from the DIC mass balance to daily rates to evaluate the correlation with daily rainfall and the extent groundwater flux could predict stream pH. Hourly fluxes were aggregated to daily rates as the sum of hourly fluxes, daily rainfall the sum of hourly rainfall, and pH as mean daily pH. To evaluate univariate control of CO2 fluxes by discharge, we used the Kendall rank correlation, τ, on daily aggregated data. We evaluated the relationship of groundwater inputs of CO2 as a driver of stream acidification events in the Taconazo. All analysis, statistics, and figures were made in R v4.0.3 (R Core Team 2020).