Regional-scale lineaments were mapped throughout the Okanagan Basin using detailed aerial (ortho) photos and LANDSAT7 Thematic Mapper multispectral panchromatic imagery (near-infrared band 4). Additional details on how the data were processed are provided in Voeckler and Allen (2012).

Airborne LiDAR digital elevation data were acquired in August 2016 (zero snow cover) using a Reigl Q1560 scanner. The average point density for last returns was 10.33 m⁻², with a horizontal accuracy of 0.3 m and a vertical accuracy of 0.15 m, both reported at a 95% confidence interval. The DEM was delivered as a GeoTIFF file at a 1-m resolution.

A suite of borehole geophysical logs, including capacitive resistivity and normal resistivity, single point resistance, magnetic susceptibility, temperature, and full wave form sonic (tube wave amplitude and variable density), was acquired over a 2-day period in June 2009. Two logging runs were acquired; a down run (logging while the probe is going down the drillhole) and an up run (logging while the probe is coming up the drillhole). This procedure provided a means of evaluating the data quality and repeatability, while also acting as a check on any drift characteristics of the sensors. Measurements were acquired as frequencies and were converted into their respective quantitative units during subsequent processing. Selected logging results for W1 and W2 are discussed in Section Evidence of Fracture Zones From Lineament Mapping and LiDAR Data.

Short duration, constant discharge pumping tests were conducted at W1 and W3; step tests were done prior to determine optimum pumping rates for the constant discharge tests and are not discussed further herein. W1 was pumped at a constant rate of 0.06 L/s for 8 h (480 min). Then, the pump was turned off and the recovery response was monitored for 30 min (achieving ~ 90% recovery). Water levels were measured both in the pumping well, W1, and in the adjacent well observation well, W2. W3 was pumped at a slightly higher pumping rate of 0.08 L/s for ~8 h (480 min), with a 30-min recovery period. As there is no other well close to W3, no observation well data were obtained for this test. Drawdown was measured in each well using a pressure transducer datalogger as well as manually with a water level tape.

The pumping test data were analyzed using different analytical methods (e.g., Theis, 1935; Cooper and Jacob, 1946; Gringarten et al., 1975), and recovery data were analyzed using the Theis Recovery method (Theis, 1935) in order to calculate the hydraulic properties of the aquifer. The software AquiferTest Pro (Waterloo Hydrogeologic Inc., 2021) was used for the analysis.

Water samples were collected from the groundwater wells (using a bailer at different depths), the soil piezometers, the stream at different locations, as well as from rainfall in summer 2007 and snow in winter 2008 (Figure 2). Temperature, pH, electrical conductivity (EC) and redox potential (Eh) were measured in the field. Half of each sample was filtered and acidified to a pH of 2 for cation analysis, while an un-acidified portion of sample was set aside for alkalinity titrations (average of 2–3 titrations), and the analysis of anions and δ²H and δ¹⁸O isotopes in water. Chemical analysis was done in the Geochemistry Lab in the Department of Earth Sciences at Simon Fraser University. Anion concentrations were measured using Ion Chromatography (IC) and major and minor cations using an Inductively Couple Plasma Atomic Emission Spectroscopy instrument (ICP-AES). Stable isotopes of oxygen and hydrogen in water were analyzed using a mass spectrometer at the Environmental Isotope Lab at the University of Waterloo.

MIKE SHE (Danish Hydraulic Institute (DHI), 2007) was used to simulate the hydrology of UPC 241. MIKE SHE is a fully distributed hydrologic modeling code that can simulate actual evapotranspiration (AET), overland flow, one-dimensional (1D) unsaturated flow, and three-dimensional (3D) variably saturated groundwater flow. Rivers, lakes, and other channels are simulated using the MIKE HYDRO River module (1D hydraulic model), formerly MIKE 11, which is coupled to the MIKE SHE model via link nodes. Further details on the comprehensive modeling capabilities of MIKE SHE can be found in the user manual (Danish Hydraulic Institute (DHI), 2007).

Figure 2 shows the location of four lineaments (black dashed lines) within UPC 241 that were mapped using orthophotos and Landsat imagery. These lineaments are inferred to be related to fracture zones. It is noted that UPC 241 is located in a region with relatively low lineament density within the Okanagan Basin; lineament density increases toward the central Okanagan valley (Voeckler and Allen, 2012), where the trace of the Okanagan Valley Fault Zone (a detachment fault), is under the Okanagan Lake and Okanagan Valley (Figure 1). The dips of these fracture zones cannot be determined based on the lineament data alone given that lineaments are visible as two-dimensional features on the landscape. The strike directions of the lineaments at UPC 241 are consistent with those mapped across the basin, and four different sets of lineaments were found to be statistically related to the four sets of fractures mapped in outcrop (Voeckler and Allen, 2012). Average dips for each of the four sets range from 11° to 58°.

Figure 4 shows the hillshade image for UPC 241 from the LiDAR data. Clearly evident is a network of elongated troughs (Figure 4A), which are interpreted to correspond to fracture zones (Figure 4B) where preferential erosion has taken place. Several of these troughs were identified in the lineament mapping (blue lines in Figure 4B), while two additional fracture zones in proximity to W3 (pink lines in Figure 4B), which were not detected by the lineament mapping, were subsequently identified in the hillshade image.

Figure 5 shows selected geophysical logs for W1 and W3: magnetic susceptibility, capacitive resistivity, tube wave amplitude and full waveform log presented as variable density log (VDL). The magnetic susceptibility log is primarily used for lithology identification, but it can also be used to map alteration zones where magnetic minerals have been altered to non-magnetic minerals. If fluid flow occurs in porous or fractured rocks, it offers an oxidizing environment which may alter magnetic minerals to non-magnetic minerals and, therefore, show as lower magnetic susceptibility zones. Lower resistivity within a rock formation is often an indicator of a discrete fracture or fracture zone, since fracture are more porous (if open) and hence exhibit low resistivity (if water-filled). Variations in resistivity in crystalline rocks are primarily a function of pore water content and salinity. Low tube wave amplitudes are often exhibited by porous fractured rocks given their low density compared to unfractured rock.

In W1, the zones of low resistivity, especially the ones below 30 m correlate well with the tube wave amplitude log and are indicated as low amplitude blue color zones in the VDL on the right of Figure 5A. The horizontal purple bars (second last column) indicate the fracture zones whose characteristics are depicted in the resistivity, susceptibility, and tube wave amplitude logs. At around 28 m, the resistivity is very low, but it is not indicated on the tube wave amplitude. The decrease of the tube wave amplitude in crystalline rocks has been correlated to permeable fractures, which suggests that the fracture zone around 28 m is not permeable. Note that the logs for W2 are essentially identical to those for W1 (not shown) given the close proximity of these wells (~3 m separation).

In W3, most of all the lower resistivity zones (Figure 5B) correlate very well with the tube wave amplitude log and are indicated as low amplitude, blue color zones in the VDL. Several fractures are identified (see fracture zone column). At shallow depth around 6 m, there is a large fracture zone about 4 m wide. In general, the fracture zones in W3 seem to have lower tube wave amplitudes than those in W1 (and W2), suggesting that the bedrock around W3 is more permeable. The differences in fracturing between the wells likely accounts for the variations in well yield. W3 has a relatively high yield at 0.761 L/s compared to W1 and W2 (0.13 and 0.19 L/s, respectively). If fracture intensity, as observed in these borehole logs, is assumed to be associated with proximity to a major fracture zone, then the low flow rates and the few fractures in W1 and W2 suggest that neither well is close to a major fracture zone, while W3, which appears to have wider and more frequent fracture zones, may be located closer to a major fracture zone.

Further evidence of the influence of significant fracturing at W3 compared to W1 and W2 is the borehole temperature logs (Figure 6). The temperature-depth profiles for W1 and W2 are virtually identical due to their close proximity, and the temperatures are virtually isothermal between 23 and 30 m, suggesting high vertical fluid flow rates in this zone, as illustrated by the orange inflow and outflow arrows on the graph. Temperature decreases slowly below 33 m depth. It is noted that fluid conductivity was relatively constant (20 μS/cm) above a depth of 32 m, but increased steadily with depth below 32 m, reaching 120 μS/cm at approximately 40 m depth. In W3, temperatures are significantly warmer than those in W1 and W2. W3 was flowing at the time of logging. This suggests there is at least one fracture zone at depth that has sufficient head to generate the flowing conditions, as represented by the green vertical arrow in the graph. However, there may be additional entry and exit intervals along this borehole.

Figure 7 shows the log-log graphs of the head drawdown vs. time and the first derivative of drawdown with time for the two constant discharge tests conducted at UPC 241: one test at W1 with observations for W1 (Figure 7A) and W2 (Figure 7B), and the other test at W3 with observations for W3 (Figure 7C). Derivative plots, commonly referred to as diagnostic plots, were used to distinguish between different response types: 1) radial flow (flow is horizontal and radial toward the well) and 2) linear (flow is one dimensional and linear in a vertical plane toward the well) (Renard et al., 2009). Derivative plots can also be used to identify spherical flow and borehole storage.

A uniformly fractured aquifer often has a radial flow response that can be represented by a classic This curve on a log-log plot (shown as the blue curves in Figure 7 for comparison) and the derivative of drawdown is constant (horizontal line) under radial flow conditions. A linear response is indicated by a more rapid initial response compared to Theis and the drawdown and its derivative are represented with straight lines on a log-log plot (Gringarten et al., 1975). Linear flow commonly occurs due to the presence of a vertical to sub-vertical discrete fracture (Gringarten et al., 1975) or wider fracture zone or dyke (Boehmer and Boonstra, 1987) that intersects, or is in close proximity to, the well (Allen and Michel, 1998). Borehole storage (BHS), which often occurs during pumping in low yielding rock units, is commonly identified on a log-log plot for the pumping well as a hump in the derivative plot at early time, and the drawdown curve has a slope of 1.

The analysis of pumping test data requires a conceptual model be identified based on the subsurface geological and the boundary conditions so that an appropriate analytical method can be selected for the analysis of the test data. The conceptual model for this study includes a network of vertical to sub-vertical fracture zones passing through a fractured bedrock unit. Based on the lineament mapping/LiDAR data and the borehole geophysics, one of these fracture zones is thought to pass close to W3. Accordingly, the possibility of linear flow influencing the pumping test at W3 was anticipated. Nevertheless, there may be other explanations for the responses to pumping at W1 and W3. The following is one interpretation. Other possibilities are explored later.

For the test at W1, drawdown data suggest a short period of BHS in the first min of the test (Figure 7A). During this time water is removed mostly from the borehole and not from the formation. Flow appears to be linear from 1 min to ~20 min. From 20 min to the end of the test at 480 min, the response appears radial (approximately a horizontal line in the derivative plot).

Figure 7B shows the test data for the observation well W2 during the pumping test in W1. Since W2 was not pumped, there is no BHS. The response is delayed slightly (by ~ 1 min) due to the 3 m separation of the wells. W2 also shows a brief period of linear flow from 1 to 20 min. Compared to W1, the radial flow period is much better defined in W2; the derivative curve is nearly horizontal.

Figure 7C shows the test data for the pumping test in W3. The drawdown curve and its derivative suggest BHS at the beginning of the test up to about one min. After that, linear flow dominates until the end of the test at 480 min.

The classical radial flow methods, Theis (1935) and Cooper and Jacob, 1946, were used to analyze the test data over the radial flow period in W1 and W2, and for the later portion of the test data in W3 despite radial flow not being observed in that test. Estimates of the hydraulic properties from both methods are shown in Supplementary Table 5. Note that storativity cannot be estimated for the pumping well. Hydraulic conductivity (K) is estimated by dividing transmissivity (T) by the representative aquifer thickness (here, the open hole interval is from the base of the well casing to the bottom of the borehole). Due to the interpreted linear response at W3, a linear flow model for pumping wells (Gringarten et al., 1975) was also used to analyze the pumping test data for W3. The entire curve was used in this analysis because it fit the data very well. The estimated properties from the Gringarten et al. methods represent the aquifer, not the fracture zone, so they can be compared with the Theis and Cooper-Jacob estimates.

The hydraulic properties from the pumping test results are very similar for W1 and W2, regardless of the analytical method used to analyze the data (Supplementary Table 5). K values range from 1.1 × 10⁻⁷ to 1.4 × 10⁻⁷ m/s. The K values calculated for W3 using the radial flow models range from 9.8 × 10⁻⁷ to 2.2 × 10⁻⁶ m/s. The Gringarten et al. (1975) method gave a K value of 1.1 × 10⁻⁶ m/s for W3. The recovery graphs are not shown, but the estimates of the hydraulic properties using the Theis Recovery method can be found in Supplementary Table 5. The results for the recovery tests are generally consistent, although the K values are very slightly higher in W1 and W2, and intermediate in W3 compared to the values from the pumping tests. Across all tests, the geometric mean K for W1 and W2 is 1.3 × 10⁻⁷ m/s and for W3 is 9.8 × 10⁻⁶ m/s.

The interpretation of pumping test data is often ambiguous because multiple conceptual models can give the same response. Uncertainties in our interpretation include 1) whether BHS occurs in the first min at both W1 and W3; 2) whether the linear response observed in W1 and W2 is due to a vertical/sub-vertical fracture intersecting those wells, or whether the response was influenced by the presence of W2 itself; and 3) whether the linear response at W3 is due to a nearby fracture zone or whether aquifer heterogeneity (e.g. a gradual decrease in K away from W3) causes the linear response. Ultimately, uncertainty in the interpretation of the pumping test data does not allow us to confirm that there are vertical / sub-vertical fractures influencing the tests; however, the high fracture intensity observed in the geophysical logs for W3 and the proximity of lineaments suggest that the conceptual model is reasonable.

Supplementary Table 6 provides the water chemistry analysis results. Water chemistry parameters vary within a relatively small range as might be expected in this headwater catchment. The pH (2008 values) ranges from 5.8–7.3 (ave. 6.5). The redox potential (Eh) ranges from 337–602 mV (ave. 445 mV), even in the deep groundwater wells, suggesting contact with the atmosphere. The low electrical conductivity (EC) ranges from 34–144 μS/cm (ave. 59 μS/cm), reflecting the low total dissolved solids (TDS) of the waters, which ranges from 27–106 mg/L (ave. 51 mg/L).

Figure 8 shows a Piper plot of select water chemistry data for 2008, as well as for the bedrock wells in 2007. Most waters have a Ca-Mg-HCO₃ or Ca-Na-HCO₃ water type. Almost all samples plot in a close cluster on the Piper plot, except a few. For example, snow samples from UP9 and UP10 located at high elevation in the catchment have elevated SO₄ compared to CP-1 at low elevation. Although the TDS values are low for all snow samples so small differences in concentrations can be accentuated. P3 has elevated Na relative to the other soil piezometers. While P3 is among the deepest piezometers (~ 100 cm), it is not substantially deeper than the other ones (depths range from 55–119.5 cm). Piezometer P11, sampled on three consecutive days in summer 2008 (P11-1,−2,−3), showed increasing relative Na concentrations. Bedrock wells W1 and W2 sampled in summer 2007 shortly after drilling have slightly higher relative Na concentrations than the same wells sampled in 2008. Samples from W3 (depths of 21 m and 30 m) have a higher relative Na concentration as well as a higher EC (Ca and total alkalinity are the highest of all samples) compared to all wells and compared to the same well in 2008. Despite the varied sampling locations within the stream network (SW1-6), the water chemistry is very similar. Sites SW1, 2, 4 and 6 have a Ca-Na-HCO3 water type, while SW3 and SW5 have a Na-Ca-HCO3 water type.

Figure 9 shows a plot of δ¹⁸O vs. δ²H for samples collected in 2007 and 2008. Rain and snow samples were used to construct a local meteoric water line (LMWL). Three rain samples (shown on the inset plot in Figure 9) have less depleted values ranging from −6 to −10 ⁰/₀₀ for δ¹⁸O and −62 to −86 ⁰/₀₀ for δ²H (Supplementary Table 6), reflecting the higher temperature during summer. Snow samples have more depleted values (shown on both the main plot and the inset plot in Figure 9). In summer 2007, samples from W1 at three depths (12 m, 27 m and 46 m) plot close to each other, slightly above the LMWL, the single sample from W2 (30 m) is more enriched (and possibly affected by evaporation), and the samples from W3 at 30 m and W3 at 21 m have a more depleted isotopic composition. In summer 2008, W1 and W2 have similar isotopic compositions, while the sample collected from W3 when it was flowing artesian in 2008 is more depleted. There is some minor deviation of the isotopic composition of all three wells between 2007 and 2008. W1 and W2 were also sampled in winter 2008 (data not shown) and plot among the main cluster of samples on the graph, so, there does not appear to be a strong seasonal effect on the isotopic composition of the bedrock wells.

The stream samples (SW1 to SW8) plot mostly along the LMWL, with SW4, SW5 and SW8 having slightly more depleted values despite their location at lower elevation in the catchment. SW5 is from the adjacent catchment UPC240 (see Figure 2). Soil piezometers P3, P4, P5 and P10 sampled in 2007 plot as a cluster above the LMWL with less depleted isotope concentrations compared to P1, P2 and P7. The variation in the isotopic composition does not appear to reflect logged vs. unlogged areas, elevation difference, or piezometer depth. Soil piezometers P3, P4 and P11 sampled in 2008 plot in the main cluster, while P12 plots lower on the line and P10 higher on the line. Interestingly, P12 is located nearby W3 in the valley at low elevation, and P10 is at moderate elevation (Figure 2).

Figure 10 compares the observed and simulated groundwater levels in W2 and W3 over the full simulation period (October 1, 1990 to July 1, 2019). As mentioned in previously, W2 continues to be monitored, while data are only available for W3 from 2007 to 2010. Considering first the results for W2 (Figure 10A), the model spin up period is very different among the four models: Model A (1999), Model B (2009), Model C (1992), and Model C (1992). Models A, C and D are dynamically stable for a much longer period compared to Model B, which took a long time to spin up and arguably may not be dynamically stable even at the end of the simulation. The long spin up period for Model B is discussed later. Comparing the model fits for W2 near the end of the simulation (from 2014–2019), the simulated seasonal amplitude of the groundwater level is slightly overestimated in Models A, C and D. The average observed amplitude is approximately 7 m, while simulated amplitudes are slightly lower. The timing of peak and low groundwater levels is rather poorly simulated by all models. Observed groundwater levels rise rapidly in October due to fall rains, continue to rise gradually until peaking in June, and then decline over the summer, while the simulated groundwater levels begin rising in late July and peak in January. All of the models produce a visually good fit for W2, although less so for Model B.

At W3 (Figure 10B), the simulated groundwater levels attain dynamic equilibrium very quickly, within about 5 years. This is because W3 is located at low elevation in the catchment where the water table is relatively shallow year-round. The four models, however, show very different average groundwater levels. Average annual groundwater levels are overestimated slightly in Model C, slightly underestimated in Model D, and greatly underestimated in Models A and B. The overall simulated seasonal amplitudes are similar among the four models, with slightly higher amplitudes compared to observed in Models A and B, and nearly identical amplitudes in Models C and D. The timing of the groundwater level responses is very close to observed for all models; however, the overall shape of the response in Models A and B is very peaked compared to the smoother observed response. Overall, Models C and D are much better at reproducing the groundwater level response in W3 than Models A and B, suggesting that the inclusion of the fractured zone and/or a weathered zone is important.

The model fit statistics for the period October 7, 2004 to October 6, 2019 for all observed data (W1, W2, W3, P1-P15, 4 snow stations, and streamflow) are provided in Supplementary Table 7. Four measures of fit are provided: root mean squared error (RMSE), normalized root mean squared error (NRMSE), R correlation, and Nash-Sutcliffe R², equivalently, the Nash-Sutcliffe Efficiency (NSE). The statistical fits for the various piezometers are mixed, varying from model to model for different piezometers. Total snow storage is very well fit at all four snow stations as is streamflow. Streamflow is overestimated by all models compared to observed streamflow by approximately 4% on an average annual basis. However, the model overestimation is primarily associated with peak flow overestimation in the months of April to June (or July). For the remainder of the months, streamflow is underestimated (as illustrated for Model D in Supplementary Figure 2). Overall, Models C and D are better at predicting the observed data than Models A and B.

Figure 11 compares the simulated depth to the top of the saturated zone (or water table) on June 18, 2010. Streamflow peaks in May or June at UPC 241, thus, June 18 represents a date when groundwater recharge would be high. The year 2010 was chosen because it is roughly midway between the end of the model spin up period and the end of the simulation. Positive numbers indicate that the water table lies above the top of the bedrock (i.e. within the soil zone), and negative numbers indicate that the water table lies at some depth in the bedrock. In all models, the depth to the water table is shallower along the streams at mid- to low elevation and surrounding the permeable vertical fractures. Presumably, because the fracture zones are closely associated with the stream network, the fractures zones may be important for maintaining baseflow in the streams during the summer. However, the water chemistry and stable isotope results for W3 are distinct from the stream samples, suggesting that a primary function of the fracture zones is the conveyance of water (i.e. snowmelt) from high elevation areas to lower elevation areas.

The daily water balance values from MIKE SHE were summarized and output every 30 days, and then averaged annually according to water year (WY), which starts October 1 and ends September 30. The annual water balance results over a 5-year period (2005–06 to 2009–10) for all models are shown in Table 1. Only the main water balance components are reported; changes in canopy, snow, overland, unsaturated and saturated zone storage, river to groundwater, and error are excluded as the average values are generally low (< ~ 1%) compared to the other components, although they do account for discrepancies on an annual basis. Consequently, the water balance items do not necessarily total the precipitation. Precipitation (P) varied annually, from 690 mm in WY 2006-07 to 802 mm in WY 2005-06. Similarly, evapotranspiration (ET) varied interannually between the four models, but not substantially so, with average values ranging from 338 mm (45% of P) to 364 mm (49% of P). Overland flow (OL) flow to rivers also varied interannually and between models, with average values ranging from 299 mm (40% of P) to 363 mm (49% of P).

In the water balance, saturated zone (SZ) to river transfer, which is the groundwater contribution to baseflow, showed almost no interannual variability, but there was a significant difference between Models A and B with a SZ to river transfer of ~1% of P compared to Models C and D with 6–7% of P; Models C and D are the two with the weathered zone. In Models A and B, the simulated baseflow is ~1% of simulated streamflow, while in Models C and D, the baseflow is 14%. Thus, the weathered zone plays an important role in maintaining baseflow in the streams. Overall, the only significant difference between models was the groundwater contribution to baseflow.

Recharge is defined as the transfer of water from the unsaturated zone (UZ) to the SZ. Recharge is not included as part of the total water balance, but rather is included in the detailed water balance results for each of the UZ (reported as a loss of water) and the SZ (reported as a gain of water). Transfers occur daily and are summed for every cell of the watershed. So, if groundwater seeps at a cell where the water table intersects the ground surface, and then water flows overland, only to enter the unsaturated zone at another cell, this would be counted as recharge. Consequently, the annual recharge numbers are misleading. Shown in Table 1 is the net annual recharge (the sum of the positive and negative numbers); the average recharge for all water years ranges from 304 mm (Model C) to 402 mm (Model A). Clearly, the numbers far exceed plausible recharge because the sum of ET, OL flow to river and SZ to river together make up 95–97% of the water balance. Thus, OL flow to river and recharge to SZ interact and cannot be clearly distinguished from each other.

Figure 12 compares monthly recharge for the year 2010 for the four models; also shown is total monthly precipitation. Positive recharge numbers indicate that water is transferred from the UZ to the SZ (true recharge), while negative numbers indicate a transfer from the SZ to the UZ. All models have consistent timing of recharge, initiating in either May or June, reaching a peak in July or August, and then declining, such that by October or November, recharge is negative. Models A and B (with no weathered zone) have high negative recharge through the late fall and winter, suggesting that there is significant transfer of water from the SZ to the UZ during this time. Once a weathered zone is introduced (Models C and D), significantly much less water is transferred back to the UZ. The abrupt change in hydraulic properties at the weathered zone – bedrock boundary, which is anywhere from 10 m to 14 m depth depending on soil thickness, likely accounts for the lower transfer. Interestingly, the weathered zone – bedrock contact depth coincides with the minimum historical groundwater level in W2 (~ 13 m). Therefore, the transition between the SZ and UZ, at least a higher elevation in the watershed, is occurring near this contact. Additionally, the high precipitation rates in November, December and January do not translate into any recharge, although water is being added to subsurface storage during these months (not shown in the annual water balance). This means that water is infiltrating the UZ and adding water to storage, but the water has not yet percolated deep enough to reach the water table; therefore, no recharge is recorded in the model.