A hydrological monitoring network was installed at Otter Park to monitor groundwater level, streamflow, and vertical flux across the streambed (Figure 3). The network includes one pumping well (PW01) and three groundwater monitoring wells (MW01, MW02, and MW03) drilled into the surficial aquifer adjacent to the site. Only MW02 was monitored in this study. It is screened in very coarse sand from 16.8–18.3 m depth. Two hydrometric stations were installed in 2019: HS2 was installed in a pool habitat in the intermittent upstream segment of the site and HS5 was installed in a downstream perennial riffle.

The hydrometric stations are maintained by the BC Ministry of Forests, Lands, Natural Resource Operations and Rural Development (FLNRORD). Discharge measurements were made by FLNRORD staff following standard procedures (see Resources Information Standards Committee, 2018) to develop station-specific rating curves to calculate streamflow from the continuous stage data. Rating curve development was done using AQUARIUS software by FLNRORD staff. Seven nested instream piezometer pairs (i.e., one deep and one shallow) were installed, with the deeper piezometers (denoted by lower numbers in Figure 3) screened at depths of around 1.05 m and the shallow piezometers (denoted by higher numbers in Figure 3) screened at depths of around 0.4 m to create a separation of around 0.65 m for each pair. Piezometer D37 was knocked down during winter high flows, so D36/D37 data are omitted from this study. Solinst Leveloggers were deployed in each of the instream piezometers, and Onset HOBO loggers were installed in the groundwater wells and hydrometric stations to log water levels and temperatures continuously. Loggers were hung from trees to measure barometric pressure. All instruments were set to record data continuously from 1 May 2020 to 16 October 2020 at 15-min intervals with some short interruptions due to data downloads.

BMI sampling was carried out to examine spatial and temporal variability in BMI communities as flows declined over the summer. BMI samples were collected at 12 sites within the Otter Park study reach (Figure 3). Six sites were assigned to riffle habitat (R1-R5) and six to pool habitat (P1-P3) to ensure an equal number of locations in each habitat type. In the upstream area where the wetted width was greater, two BMI samples were collected along each transect; the sample location to the north is denoted with an A, while the south location is denoted with a B (Figure 3). The six sites within each habitat type were evenly distributed among locations with relatively high streambed elevation and relatively low streambed elevation in case water preferentially pooled in low elevation sections as flow declined. The upstream pool habitat unit, extending from P3 downstream to P1, had lower streambed elevations than the upstream riffle habitat unit, which extended from R5 downstream to R3 (Figure 3). The downstream riffle habitat unit, extending from R2 downstream to R1 had lower streambed elevations than both upstream habitats units (Supplementary Figure 1).

BMI samples were collected at each site once per month from June to September in 2020 (June 11, July 14, August 10, and September 17). The sampling events were timed to occur roughly 30 days apart to allow for macroinvertebrates to recolonize substrates disturbed during each sampling event. Each visit began at site R3 and progressed upstream to not disturb other sites during sampling. The two furthest downstream sites, R1 and R2, were sampled last because the channel length and complex channel features between R2 and R3 likely prevented disturbance from upstream activity.

For each BMI sample, the 30.5 × 30.5 cm frame of a 500 µm mesh Surber sampler was centered next to the site marker and the substrate within the frame was disturbed by hand vigorously for 2 minutes to dislodge macroinvertebrates into the net on the downstream end. After 2 minutes, the net was removed, and its contents were placed in sample containers and preserved with isopropyl alcohol. In September, flow was too low in pools to direct dislodged macroinvertebrates into the net, so the substrates were scooped by hand into the net. Various techniques may be used to sample stagnant pools (e.g., Bonada et al., 2006; Chester and Robson 2011), but modifying the Surber procedure allowed for consistency in sample area, substrate disturbance mechanism (hand), and the steps involved in removing specimens from the net. Macroinvertebrate specimens were stored and identified within 2 weeks after their collection to the family level using a taxonomic key (Bouchard, 2004).

Habitat variables were measured after each macroinvertebrate sample to characterize spatial and temporal variability in habitat conditions. First, discharge measurements were made using the Sontek FlowTracker2 along transects crosscutting cutting each site. Depth and velocity were measured in equal increments along each transect for a total of approximately 23 measurements per transect. By keeping the number of measurements relatively constant for all transects, the spacing of measurements was smaller for transects with smaller wetted widths. Discharge for each transect was calculated in the FlowTracker2 desktop application using the mean-section method (Sontek 2019). After each transect was completed, dissolved oxygen, pH, and conductivity were measured with a YSI 600 XLM multiparameter water quality sonde at the macroinvertebrate sampling points.

The depth and velocity measurements made on either side of the Surber net were averaged to assign a two-point average to each macroinvertebrate sample. Velocity data were missing for sites R3 and R4 in September because the water was not deep enough to measure velocity with the FlowTracker2. The missing values were estimated by multiplying the measured September depths by the average velocity:depth ratios for each corresponding site. Depth measurements were not available to replace missing June data for P2A and P2B.

A substrate survey was conducted at each site in June. Turbid water at some of the upstream sites prevented visual assessment of the substrate composition. Instead, the 100-stone count used by the Canadian Aquatic Biomonitoring Network (CABIN) to characterize substrate at the mesohabitat scale was adapted to characterize substrate for smaller habitat sizes (Environment Canada 2011). Twenty random stones were selected from the 50 cm radius surrounding each site marker and the intermediate axis of the stone was measured and recorded. The intermediate axis is perpendicular to the longest axis of a stone and the measurement is taken at the point resulting in the longest intermediate axis measurement (Environment Canada 2011).

To examine potential differential temperature patterns in riffle and pool habitat, temperature recorded by the HOBO loggers installed at HS2 and HS5 were used to represent stream temperature in pools and riffles, respectively. Groundwater temperature, air temperature, and canopy coverage were taken into consideration to support the analysis of stream temperature. Groundwater temperature was represented by the temperature recorded by the HOBO logger in MW02. Canopy coverage for the area around the hydrometric stations was visually assessed and recorded as a percentage of the stream section adjacent to each station with overhead canopy cover (Environment Canada 2011). One HOBO TidbiT v2 temperature logger was attached to a tree adjacent to the site to monitor air temperature.

Data manipulation, visualization, and analysis were carried out in R Version 4.0.3 through the integrated development environment R-Studio Version 1.3.1093 (R Core Team 2020). The R packages ggplot2 (Wickham 2016), ggfortify (Tang et al., 2016), RcolorBrewer (Neuwirth 2014), dplyr (Wickham et al., 2020), tidyr (Wickham 2021) were used frequently.

Hydrological monitoring network data were analyzed to characterize variability in groundwater level, streamflow, and vertical flow. Daily precipitation data recorded at the Abbotsford International Airport climate station were compared to mean daily groundwater level, discharge at HS2 and HS5, and vertical flow at each piezometer pair. Discharge data from the two hydrometric stations, HS2 and HS5, were plotted on a pseudo-log scale using the R package scales (Wickham and Seidel 2020) to characterize flow variability, specifically in relation to low flows and zero-flow events. Hydraulic gradients across the streambed were calculated at each 15-minute time step for each piezometer pair by taking the difference in water levels between the nested pair and dividing by the difference in their depth. The vertical flow across the upstream section of streambed between transects R5 and P3 (Figure 3) was estimated for each 15-min time step from the instantaneous hydraulic gradient values calculated from instream piezometer water level data (See Supplementary Material for gradient calculations), the streambed area (510 m²), and the average streambed hydraulic conductivity (2.5 × 10⁻⁵ m/s), which was estimated from slug tests carried out in the instream piezometers during a previous phase of field investigations at the site (Allen et al., 2020).

Hydraulic habitat variables (depth, velocity, and substrate) were analyzed to characterize differences between riffles and pools, and variability in habitat conditions over the summer. For both depth and velocity, a two-way analysis of variance (ANOVA) was used to analyze the potential effects of habitat type, month, and the interaction between habitat type and month. Velocity was transformed using Tukey’s Ladder of Powers transformation through the R package rcompanion (Mangiafico 2016). One-way ANOVA was used to test for differences in substrate size among sites.

A Welch’s Two Sample t-test was used to see if substrate differed between all sites in pools and all sites in riffles. The df and p are reported in brackets for t-tests.

For all ANOVA performed in this study (including those detailed in Sections 2.7 and 2.8.1), results are provided in brackets starting with the F-statistic (F), with the degrees of freedom (df) and sample size (n) as the first and second subscripts, respectively (e.g., F_df,n). The p-values (p) are reported after F in each set of brackets (e.g., F_df,n = 1, p < 0.05). Post-hoc Tukey’s Honestly Significant Difference (HSD) tests were used to test for significant differences between all possible levels of each factor (i.e., month, site). Tukey’s HSD results are reported in brackets, with p < 0.05 indicating significant differences and p > 0.05 indicating no significant differences. Levene’s test of homogeneity of variance was performed through the R package car (Fox and Weisberg, 2019), and diagnostic plots were visually inspected to ensure that the assumptions of the tests were met (Supplementary Figures 2–8).

Water quality variables (DO, pH, EC, and temperature) were analyzed to characterize spatiotemporal variability and assess their potential influence on BMI communities. Time series of stream temperature, groundwater temperature, and air temperature were plotted. Data missing from June due to a malfunctioning DO probe were replaced with DO recorded 13 days after the original sampling event. Discrete DO and pH data, and stream temperature time series data from HS2 and HS5, were evaluated against water quality guidelines (Canadian Council of Ministers of the Environment (CCME), 1999; ECCS 2019). EC was evaluated against the typical values of EC found in freshwater streams. If some values of a water quality variable exceeded the guideline, a two-way ANOVA test was performed to see if the variable differed between sample sites or months.

A generalized additive model (GAM) was built to assess the potential effects of habitat variables on macroinvertebrate abundance. GAMs are generalized linear models in which the response variable depends on the sum of smooth functions of one or more predictor variables (Wood 2017). Potential complex non-linear relationships between predictor variables and response variables justified the use of GAMs over linear models.

Predictor variables included in the analysis were EC, pH, DO, velocity, depth, and substrate size. For depth and velocity, the two-point averages calculated for each sample were used in the GAM. Substrate was incorporated as a combined substrate index designed to better capture the variability in substrate sizes (Jowett et al., 1991). To calculate the substrate index, stones comprising each 20-stone sample were classified using the modified Wentworth particle size scale (Bovee and Milhous 1978) adjusted to include a fine gravel category (Jowett et al., 1991). The weighted percentages of each size class were then summed to get the final substrate index (Jowett et al., 1991): Substrate index =   0.08  bedrock%  +   0.07  boulder% + 0.06  cobble% + 0.05  coarse gravel%  +   0.04  fine gravel%  +   0.03  sand% . (2)

The GAM for macroinvertebrate abundance was fit using the R package mgcv. The functional form of the relationship between the predictor variables and the response variable was defined as: Y   ∼   s ( x 1 ) +   s ( x 2 ) +   s ( x 3 ) + s ( x 4 ) +   s ( x 5 ) + s ( x 6 ) (3) where Y is macroinvertebrate abundance, and   x 1 , x 2 , x 3 , x 4 x 5 , x 6 are velocity, depth, substrate index, DO, pH, and EC. BMI abundance was cube-root transformed to meet the model assumptions. The GAM was fit using restricted maximum likelihood (REML) smoothing parameter estimation, which is considered to be the most stable method (Wood 2017). The distributional family, which defines the distribution of the response conditional on the covariates, was set as Gaussian. While Poisson and quasi-Poisson distributions are often used to model count data, consequences of violating these distributions are severe, whereas Gaussian models are more robust to violations (Hoef and Boveng 2007; Knief and Forstmeier 2018). The double penalty approach was used to identify important model variables (Marra and Wood 2011). Model selection was carried out using Akaike Information Criterion adjusted for small sample sizes (AICc). After each model run, the variable with the lowest p-value was dropped and this process was repeated until the AICc score was maximized, GAM results and diagnostics were plotted using the R package gratia (Simpson 2021). See Supplementary Figure 9 for GAM diagnostics.

Total monthly precipitation decreased from 111 mm in May to 44 mm in August, with July and September totals just exceeding 50 mm (Figure 4A). The frequency of precipitation events also decreased over the summer, with just a few events punctuating extended dry periods. In response to decreasing precipitation, the groundwater level in MW02 decreased by 0.4 m from 2.9 m below ground surface at the start of May to a low of 3.3 m below ground surface by mid-September (Figure 4A).

Vertical flow across the streambed, averaged over the stream segment between HS2 and R5, was highly variable over the summer (Figure 4B). Precipitation events caused abrupt increases in the magnitude of vertical flow in the downward direction (from stream to aquifer). However, the mean vertical flow was 5.7 × 10⁻⁵ m³/s and was generally in the upward direction when it was not raining. Vertical flow peaked above 2.5 × 10⁻⁴ m³/s in early September when precipitation was low and infrequent.

Over the study period, mean streamflow values at HS2 (Figure 4C) and HS5 (Figure 4D) were 0.15 m³/s and 0.015 m³/s, respectively. Streamflow at both stations generally remained above 0.1 m³/s from May to early July and then decreased to values near 0.01 m³/s by August. Flow persisted at both stations from August to September. However, flow ceased at HS2 on September 6. Flow cessation was punctuated by a few hours of measurable flow, which coincided with the macroinvertebrate sampling event carried out on September 17. Otherwise, the period of flow cessation generally extended until September 23.

During this period of flow cessation at HS2, flow persisted at HS5 and there was an increase in the magnitude of vertical flow in the upward direction in the upstream section between transect R5 and transect P3. The increase in upward vertical flow indicates an increase in the groundwater discharge to the stream. Since there was no precipitation during this period, groundwater discharge to the stream likely maintained streamflow downstream at HS5. In turn, the role of groundwater in maintaining instream flow conditions likely became more important over time as precipitation decreased.

Variability in streamflow and aquifer-stream exchange drives expansion and contraction of IRES networks and therefore dictates the timing, magnitude, and duration of drying, which controls stream biodiversity at large scales (Boulton et al., 2017). High variability in flow systems across IRES, as suggested by their large global extent, stresses the importance of research on IRES flow systems (Datry et al., 2014). Like other IRES (Gómez-Gener et al., 2021), flow paths likely differed between the flowing phase (Figure 11A) and the period of flow cessation (Figure 11B). During the flowing phase, maintenance of streamflow by precipitation and groundwater input is common in IRES (Woessner 2020). At Otter Park, the flowing phase lasted from May to early September. The average vertical flow across the streambed between transect R5 and transect P3 was frequently upward, indicating the stream was gaining groundwater. The hypothetical flow paths during this phase as illustrated in Figure 11A would resemble flow paths of a perennial stream (Winter et al., 1998; Sophocleous 2002). Shallow hyporheic flow is driven by streambed topography and interacts with local groundwater flow paths (Harvey and Bencala 1993; Sophocleous 2002).

During the period of flow cessation, the hypothesized flow system illustrated in Figure 11B was formulated by referencing the hydrometric and vertical flow data against conceptual models found in the literature. Firstly, the maintained connectivity between the upstream pool and the downstream riffle is consistent with previous conceptual models of flow systems occurring in isolated pool networks in IRES during flow cessation (Datry et al., 2014; Rau et al., 2017; Gómez-Gener et al., 2021). The persistence of flow at HS5 through the entire period of flow cessation indicates that water was at least being sourced from draining of upstream pools through the hyporheic zone. Without some source of water to replenish what was flowing out of the downstream pool and riffle at HS5, the downstream segment would have ceased to flow and likely dried. In addition to hyporheic flow from the upstream pool to the downstream pool, groundwater discharge along the stream, as evidenced by the large positive vertical flow observed for the upstream pool, also likely helped maintain flow at HS5. As the influence of groundwater is typically diluted at higher flows (Gómez-Gener et al., 2021), groundwater contribution to the stream became more important for maintaining wetted instream conditions as precipitation decreased and streamflow declined. Perennial reaches are often assumed to be valuable refuges where the persistence of ideal habitat conditions yields more diverse and rich communities (Burk and Kenney 2013; Datry et al., 2014).

Physical and Chemical habitat characteristics differed between riffles and pools. Depth decreased more rapidly in riffles than in the upstream pool. Faster drying and habitat loss in riffles is a key characteristic of habitat variability in IRES (Gomez Gener et al., 2021). Higher velocity and lower depth for most riffle sites across most months relative to pools are typical of riffle-pool differences in all stream types since depth and velocity are two variables used to distinguish between habitat types (Brown and Brussock, 1991; Wohl et al., 1993). Finer substrate is often a distinguishing characteristic of pools. But fine material may be scoured from the streambed during periods of higher flow and deposited during periods of lower flow. Since streamflow was at least an order of magnitude higher in June relative to much of the remainder of the study period, spatial heterogeneity in substrate size may have been greater in other months.

In addition to depth and velocity in riffles being more sensitive to streamflow decline relative to pools, stream temperature in the riffle was more sensitive to changes in air temperature. Maintenance of lower temperatures relative to the riffle through the period of flow cessation can likely be attributed to pool depth remaining near 30 cm (Poole and Berman 2001; Caissie 2006). The decline in both riffle and pool temperature from August to October, including over the period of flow cessation, indicates that temperature was not the only driver of dissolved oxygen decreases over the summer. Instead, low DO in September was likely the result of oxidation of organic matter in isolated pools over the period of flow cessation. Similar DO in riffles and pools is common for low gradient streams like Bertrand Creek because re-aeration by turbulent flow in riffles generally declines with channel slope (Rosenfeld 2017).

Differences in habitat conditions between individual sites suggest that habitat heterogeneity occurs at the site scale in addition to the scale of habitat type. For example, R3 and R4 were consistently higher in velocity and lower in depth than all other sites, DO was likely highest at R3 and R4 in September because their high velocity and low depth allowed for local re-aeration by turbulent flow. While turbulent flow did not appear to be an important control on DO, riffle aeration might have a greater influence on local DO concentration later in the summer when DO is already low.

In this study, BMI community composition, richness, abundance, and diversity differed between riffles and pools, but also between the downstream perennial riffle and the upstream intermittent riffle. The greater abundance of most taxa in riffles but more even distribution of chironomids between riffles and pools are not uncommon (Brown and Brussock, 1991). Baetidae and Simuliidae are significant identifiers of riffle communities because they are rheophilic (flow-loving) taxa (Scullion et al., 1982; Boulton 2003; Bonada et al., 2006; Stubbington et al., 2009). However, 87% of Baetidae and 96% Simuliidae individuals were found at the two upstream riffle sites, indicating heterogeneity at the site scale may be more important than heterogeneity between habitat types for predicting the spatial variability in the abundance of these taxa. While some baetid mayflies and black flies were found at the downstream perennial riffle sites in some months, they had a similar dominance of chironomids and gammarids as the intermittent pool sites. Consequently, richness, abundance, and diversity did not differ between the upstream pool and the perennial riffle, but all metrics were significantly higher for the upstream riffle, which contrasts the common finding that community compositions at intermittent sites are subsets of those found at perennial sites (Feminella, 1996; Price et al., 2003; Bonada et al., 2006; Datry et al., 2014). Even in September, 1 day after flow resumed following 10 days of flow cessation, richness, diversity, and abundance remained elevated at the two intermittent riffle sites. In Mediterranean climate streams in California, Bonada et al. (2006), also found that community compositions were similar in perennial riffles and pools during the flowing phase; however, this was attributed to the direct connection between the riffle and the pool. At Otter Park, the pool and perennial riffle are separated by over 100 m. Bonada et al. (2006) asserts that sites in different habitat types in the same stream reach or the same stream are more likely to be similar than communities in the same habitat type but different streams because larger scale catchment characteristics act as the first filter, while differences at the scale of habitat type act as the second filter. Therefore, larger scale filters could be responsible for the similar communities found in the upstream pool and the downstream perennial riffle. However, the consistently higher diversity, abundance, and richness at the two intermittent riffle sites still suggests high importance of habitat heterogeneity at the site scale.

The consistently low richness and diversity of BMI communities across the site likely reflects how recurring patterns of extreme variability in physical and chemical habitat characteristics associated with streamflow decline, flow cessation, and drying have sculpted the BMI community in that only taxa with optimal traits for persistence inhabit reaches near Otter Park. Each taxon present in reasonable abundance from June to September possesses traits that promote their survival through flow cessation and drying. For example, taxa present immediately following flow cessation included Chironomidae, Tabanidae, and Gammaridae, which are common inhabitants of IRES (Nimmo et al., 2015; Datry et al., 2017). During a period of flow cessation in Sycamore Creek, Arizona, Tabanidae were commonly found taking refuge in dry substratum and also took refuge below dried leaf litter (Boulton 2003). The refuges were also occupied by chironomids, which are also commonly found taking refuge in the hyporheic zone (Boulton 2003; Datry et al., 2017). Gammarids also often enter the hyporheic zone during flow cessation and drying and are early colonists when flow resumes (Vander Vorste et al., 2016; Stubbington et al., 2018). Gammarids also have strong swimming abilities, which allow them to swim upstream from reaches where flow persists to recolonize the upstream riffle when flow resumes (Poff et al., 1997). At Otter Park, these taxa likely still dominate at the downstream riffle because low flows, low dissolved oxygen, and high temperature variability are still fundamental aspects of the natural variability of habitat conditions at the downstream riffle.

Although flow cessation and drying filter out intolerant taxa across the entire reach, the high abundances of black flies and baetid mayflies suggest that the habitat preferences of certain taxa can still be important controls on spatial and temporal variability in abundance. In Canada, the abundances of baetid mayflies are strongly correlated with velocity (Monk et al., 2017). Black fly abundances have also been shown to increase with velocity (Rivers-Moore et al., 2006). Burgazzi et al. (2020) found that baetid mayflies and chironomids occupied a wide range of habitat conditions, but abundances were driven by habitat preferences. However, gammarid abundances show a negative relationship with velocity (Monk et al., 2017).

Dissolved oxygen may have contributed to higher abundances at the intermittent riffle sites following flow resumption, but our confidence in the effect of DO is low because the DO was not consistently higher than other sites like velocity was. The GAM analysis and the NMDS environmental vector overlay analysis indicated DO and velocity were likely important controls on total abundance and BMI community composition, respectively. However, the results of these two analyses should be interpreted as an indication that velocity and DO may be important drivers of BMI abundance and community composition until they have been replicated. Other factors, such as predation, food availability, habitat availability, and potential seasonal effects related to the life cycles of different taxa could influence the accuracy of the results. Regardless, four consecutive months of higher abundance, diversity, and richness at two intermittent riffle sites relative to an intermittent pool and perennial riffle likely indicate that site-scale habitat heterogeneity influences spatial variability of BMI communities even though flow cessation and drying act as a primary filter, which limits the overall taxonomic richness across the site (Bonada et al., 2006).

Although groundwater discharge to the stream likely helped maintain streamflow in the downstream riffle during the period of flow cessation, BMI communities did not increase in diversity, richness, and abundance as predicted in other studies (Datry et al., 2014). However, these findings likely highlight the immense complexity of groundwater dependence in IRES rather than a lack of ecological dependence on groundwater. For example, upward vertical flow across the streambed in the upstream pool suggests that groundwater helped maintain water level in the pool, even during flow cessation. Unravelling the complex dynamics of aquifer-stream exchanges in IRES will likely continue to reveal how changes in hydrological connectivity help organize aquatic communities and their physical habitat templates.