Streamflow and weather conditions of seven small coastal watersheds, British Columbia, Canada, 2013-2019

Hydrometeorological observations of small watersheds of the northeast Pacific coastal temperate 10 rainforest (NPCTR) of North America are important to understand land to ocean ecological connections and to provide the scientific basis for regional environmental management decisions. The Hakai Institute operates a densely networked and long-term hydrometeorological monitoring observatory, that fills a spatial data gap in the remote and sparsely gauged outer coast of the NPCTR. Here we present the first five water years (October 2013–October 2019) of hourly streamflow and weather data from seven small (< 13 km), coastal watersheds. Average yearly rainfall was 15 3267 mm, resulting in 2317 mm of runoff and 0.1087 km of freshwater exports from all seven watersheds per year. However, rainfall and runoff were highly variable depending on location and elevation. The seven watersheds have rainfall-dominated (pluvial) streamflow regimes, streamflow responses are rapid and most water exports are driven by high-intensity fall and winter storm events. Measuring rainfall and streamflow in remote and topographically complex rainforest environments is challenging, hence advanced and novel automated measurement methods were 20 used. These methods, specifically for stream flow measurement allowed us to quantify uncertainty and identify key sources of error, which varied by gauging location. Links to the complete dataset, watershed delineations with metrics, and calculation scripts can be found in Sect. 6 and 7.


70
(Environment Canada stations with gross drainage area < 25 km 2 and USGS stations with 24 h average summer flow < 1.5 m 3 s -1 ) are shown to indicate the lack of gauging stations measuring outer coast, small watersheds. Terrain relief map was created by Arriola and Holmes (2017). Regional rainforest cover is derived from the original rainforest distribution mapping of Ecotrust (1995), which reflects the rainforest subzones of Alaback (1996).  495 m a.s.l. on central Hecate Island to relatively low gradient hummocky terrain in the west. The landscape is characterized by extensive wetlands, bog forests, and bog woodlands (Green, 2014;Thompson et al., 2016), with shallow (typically < 1 m) but organic-rich soils (Oliver et al., 2017), on faulted and folded intrusive igneous bedrock 90 (primarily quartz diorite) (Roddick, 1996). Forest stands are generally short with open canopies and tree composition is dominated by western redcedar, yellow cedar, shore pine, and western hemlock (Thompson et al., 2016;Hoffman et al., 2021). Understory and wetland plants include bryophytes, salal, deer fern, sphagnum mosses and sedges.
Calvert and Hecate Islands are located within the Hakai Lúxvbálís Conservancy (protected area) and the unceded 95 territories of the Haíłzaqv and Wuikinuxv Nations. The study watersheds have no damns or diversions, no active roads, and little evidence of historic logging except in some shoreline areas which are now well forested. Fire is typically infrequent in the coastal temperate rainforest yet recent work has revealed legacies of long term cultural burning near village sites in this area . The KWO design captures the seven largest (3.0 to 12.8 km 2 ) watersheds draining into Kwakshua and Meay Channel, signified as 626, 708, 703 and 693 (Calvert Island) and 1015, 819, and 844 (Hecate Island). The watersheds are on average 68 % forested, 24 % non-forested but vegetated (including open wetlands), 4 % nonvegetated (including exposed bedrock), and 4 % covered by lakes (Table 1 ). A total of 159 lakes (> 0.02 ha) are 105 unequally distributed across watersheds (range of 6 -53 lakes per watershed). Maximum watershed elevations range from 160 to 1012 m a.s.l. and slopes are on average 31 % (Table 1).
Catchment characteristics vary among the seven gauged watersheds (Table 1 7.5% respectively). Watershed 626 is typified by multiple small lakes and a relatively large area of exposed bedrock and non-forested vegetated land (e.g. bogs, 59%). Watershed 708 is typified by one centrally located large (30 ha) 120 lake. Watershed 1015 (3.33 km 2 ) has a subdued topography in the west, characterized by a large lake (28 ha.) near the outlet, which is flanked by a steeply sloped ridge in the east (max elevation 432 m a.s.l.). The watershed has a relatively high forest cover (80%).
Each watershed has a hydrometric station near the outlet and a meteorological station at low to mid elevation. In addition, three meteorological stations have been installed at higher elevations between watershed boundaries ( Fig.   125 1, Table 2). Photos of meteorological and hydrometric stations are given in Fig. A1 and A2.

Instrumentation and data collection
Seven hydrometric and fourteen meteorological stations were installed in a tight network spanning most of the elevation gradient Fig. 1), as the study area's complex topography leads to rapid changes in 130 climatic parameters over short distances (Beniston, 2006). Station details (installation date, location, and elevation) are provided in Table 2, and sensor inventory and specifications can be found in  Air temperature, relative humidity and rain were measured at two meters above ground for all stations, except 'East Buxton' where seasonal snowpacks necessitate the precipitation gauge orifice and air temperature/relative humidity sensor to be installed at 4.5 and 4 meters above ground respectively. Wind speed and direction were measured at 10 m above the ground at 'Reference Station', 8 m at 'East Buxton' and at 5 m at all other stations. Tipping bucket rain 160 gauges (TBRG) were field calibrated semi-annually using a Field Calibration Device (FCD-653, Hydrological Services, Lake Worth, USA) with a 200 mm h -1 nozzle rate. For the field test, the FCD was emptied into the TBRG twice and the total number of tips were recorded. The current field specification required the TBRG to be within +/-6 tips of 202; where 202 tips is the expected number of tips delivered by two applications of the FCD. If the results of the field test exceeded this threshold, a correction factor was applied to the TBRG to adjust the volume per tip 165 ratio. The air temperature/relative humidity sensors were replaced every two years. The sensors were field checked twice per year in between replacements using a 'benchmark' sensor (H2SC3, Campbell Scientific, Edmonton, Canada). In addition, a second temperature sensor (CS 109, Campbell Scientific) was installed at each site as a check and back-up in case of main sensor malfunction. The measurements of the main sensors were corrected if the difference between the benchmark and the main sensor measurements exceeded the sum of their accuracies (± 0.2 170 ℃ and ± 1.6 %).

Hydrometric stations
Stream water level (stage) was continuously measured at each hydrometric station (starting fall 2013 for watershed 708 and fall 2014 for all other watersheds) and periodic discharge measurements were taken along the range of potential water levels to develop stage-discharge rating curves.

175
Station locations were selected based on channel stability and potential to measure discharge. Pressure transducers (OTT PLS-L, 2016), anchored to bedrock or boulders with cables armored with steel, were used to measure stage (0-4 m range SDI-12). The mean, maximum, minimum, and standard deviation of five second sampling intervals were recorded every five minutes. Stand-alone water level recorders (Odyssey Capacitance Water Level recorder -Data Flow Systems PTY Ltd 2016, replaced by HOBO Water Level Data Logger -U20L-04 in September 2018) were 180 installed in proximity to each pressure transducer as a backup in case of sensor or data logger malfunction. Relative location of each pressure transducer was regularly surveyed to a benchmark to help identify potential sensor movements. Stream conditions and streambed morphology were monitored through time lapse cameras pointing at each station's stream cross-section and through photos taken during maintenance visits.
Low flows, generally below 0.5 m 3 s -1 , were measured manually using the velocity-area method midsection stream flow measurements at high water levels. Therefore, a novel discharge measurement system was designed based on the salt in solution method as described by Moore (2005): salt in solution (5L water to 1kg salt) is stored in a 1000L IPC tote on site, with a pump to pre-mix the salt solution before a measurement, and another pump to transfer the solution to a stainless-steel cylinder used to calculate the volume, before it is transferred to a dumping 195 mechanism over the stream that holds up to 40L (Fig. A2). At predefined stages -based on gaps in the stagedischarge rating curve -a signal is sent to release a predetermined volume of salt solution. This volume is targeted to never exceed the most sensitive toxicity threshold of 400 mg L -1 (Moore 2004a(Moore , 2004b. Upon initiation of the salt solution pump sequence, a second command is sent to a downstream data logger to activate at a minimum two electrical conductivity sensors (Global Water instrumentation, Inc., College Station, USA; T-HRECS Fathom

200
Scientific Ltd), installed at opposing sides of the stream to capture the passing salt wave. Upon completion of the dump sequence, the EC data are transmitted via the telemetry network for data storage and processing. Recharging of the salt solution reservoir was done manually and a calibration factor (CF), needed to convert EC sensor readings to the relative salinity of the salt solution, was determined prior to refill as well as after refill, with a goal of two CFs per salt fill. To increase accuracy, a triplicate reading was taken for each CF and each salt in solution measurement 205 was matched with one set of triplicate CF values based on their sampling date. Finally, discharge was calculated using the salt volume data and relative salinity (Moore, 2005) for each EC sensor and for each CF. This yielded six discharge values per measurement, which were averaged to calculate final discharge.

Data analysis and quality control
Characterization of data quality was done by two descriptors, which were stored together with each observation: 210 data processing levels and data quality flags. Data processing levels indicate the status of data handling: unpublished, raw data are termed 'level 1', data subjected to quality control are labelled 'level 2', and 'level 3' refers to derived data products (e.g., gap-filled data). The flagging scheme consists of two tiers: the first tier includes generic flags, e.g., 'Accepted Value' (AV), 'Estimated Value' (EV), 'Suspicious Value Caution' (SVC), and 'Suspicious Value Discard' (SVD) (Henshaw & Martin, 2014). The second tier is a use-case-specific comment on 215 data quality (e.g., background events affecting data values or failed individual quality tests). Any changes or corrections applied to the data are stated in the second tier allowing data users to customize data filling specific to their research objectives. In addition to general quality checks, a quantitative uncertainty analysis of the streamflow data was performed, focusing on discharge measurement and rating curve error. The following section provides details specific to each dataset for quality assurance and quality control.

Weather data
Measured weather data includes air temperature (℃), relative humidity (%), rain (mm), total precipitation (mm), snow depth (m), wind speed (m s -1 ), wind direction (degree), and incoming solar radiation (W m -2 ). Quality assurance procedures involved a combination of visual and automated inspection of the data for inconsistencies, such as sudden large increases in measurements, as well as outliers. Table A2 details the thresholds used to identify 225 outliers, and the rate of change (ROC) used to identify inconsistencies in the data. All data, except for the snow depth data, were gap-filled using linear regression from nearby station data. For instances when air temperature https://doi.org/10.5194/essd-2021-423

230
Precipitation data required additional quality assurance procedures due to errors introduced by wind. Tipping bucket rain gauge (TBRG) data from 'WSN626', 'WSN693_703', 'WSN703_708', 'WSN819_1015', 'Hecate', 'East Buxton', 'SSN693', 'Reference Station', and total precipitation data from 'East Buxton' and 'Reference Station' were corrected for undercatch (Yang et al., 2008) using the wind speed adjustment from Legates et al. (2004). In addition to these corrections, suspect spikes in rain intensity caused by wind-induced tips were identified, initially 235 through visual inspection, but starting from February 2016 through an automated method: when three or more tips occured within a 5 s scan interval, the data were flagged assuming that 3 or more tips within 5 s would yield an extremely unlikely rainfall rate for the area (> 200 mm h -1 ). All suspect tips were removed and gap-filled. Obvious spikes in the snow data (> 1 m) that did not correspond to increases in total precipitation were removed, but recurrent gaps in the snow depth data that were caused by sensor failures were too large to be gap-filled.

Streamflow data
Stage-discharge rating curves were calculated for each watershed (Fig. 2). Rating curves were updated generally once per year, with each site having between 5 to 50 new measurements.
Stage data were quality controlled and flagged by visual inspection; any spikes or drops and changes in stage not associated with a rain event were further inspected and where necessary flagged or corrected. Time-lapse photos 245 were also used to assess the turbulence of flows affecting sensor measurements. Data gaps < 30 min were filled using linear interpolation while longer gaps were filled using linear regression of the data from the back-up standalone water level sensors that were installed in proximity of the main stage sensors. In the one exceptional case where both water level sensors were malfunctioning (Watershed 819, 1 Dec 2015-17 May 2016), the stage records of watershed 844, the neighboring watershed sharing many catchment characteristics, were used to estimate stage by 250 linear regression.
Discharge measurement uncertainties were calculated using the Interpolated Variance Estimator (Cohn et al., 2013) for the Velocity-Area method, which estimates uncertainty from the fluctuations in the flow velocity and depth profile, as well as from calibration errors in the width, depth and velocity measurements. Uncertainty estimations for the automated 'salt solution' discharge measurements were recorded in both a quantitative and a qualitative manner.

255
First, a quantitative % error was calculated by error propagation of the uncertainty in the volume of salt solution, the EC sensor resolution, and the salt solution-electrical conductivity calibration factor (Korver et al., 2018).
Measurements with uncertainties higher than 20 % were excluded from further analysis. Second, the EC data were assessed for noise and spikes, caused by turbulence at high flows; where possible spikes were removed, and noisy data were smoothed. Measurements that contained EC readings so noisy that they could not be corrected, were EC sensor site was performed by comparing calculated discharge from the two EC sensors placed within opposing sides of the river. If calculated discharge differed more than 1 %, mixing conditions were considered poor and, unless the measurement filled an essential gap on the stage-discharge rating curve, it was excluded from further 265 analysis.
Rating curves for each station were plotted using locally estimated scatterplot smoothing (LOESS) regression of the stage-discharge measurements (Fig. 2). The span widths of the LOESS fits were visually assessed and selected.
Discharge data prior to and after high-flow events were analyzed for possible shifts in the rating curves; in case of a shift, time-lapse videos and photos taken during maintenance surveys were investigated to confirm a concurrent 270 change in streambed morphology. The curves were extrapolated to minimum and maximum stage of the five-year stage time-series, estimating minimum stage to equal zero flow, and estimating maximum flow by extrapolating a power-law equation fitted on the upper section of the rating curves.
Following the methodology proposed by Coxon et al. (2015), rating curve uncertainties were quantified by plotting 95 % confidence intervals (CI) around the curves. CIs were derived from 500 curve fitting results of LOESS 275 regressions on 500 randomized sets of stage-discharge measurements and their maximum and minimum absolute measurement errors. Measurement error of a single stage-discharge set was calculated by error propagation of the discharge measurement uncertainty (described above) and the stage measurement uncertainty. The latter were calculated from the standard deviation of stage values recorded during the discharge measurement, so that measurements taken during rapidly falling or rising water levels get assigned higher uncertainties. For each stage-280 discharge measurement set on the rating curve, 500 discharge values and 500 standard deviations are predicted and combined in a Gaussian mixed model, to derive minimum and maximum absolute discharge (95 % CI). Any minimum flow extending below zero was set to zero. Calculations were done in R and the code has been made publicly available (Sect. 7). The output is a stage-discharge lookup table with values of discharge, minimum discharge (95 % CI) and maximum discharge (95% CI) for each mm of stage.  Precipitation from 2015 to 2019 averaged 3246 mm over all watersheds, however this may have been slightly underestimated as precipitation was measured as rain only, disregarding snow, at all but two locations ( Table 2).
Temperatures were generally mild with mean annual temperatures of 8.2 °C ranging from 2.6 °C in December to 14.5 °C in August. Winter was characterized by high windstorms (>10 m s -1 ) from the southeast shifting to calmer 295 winds from the south-and northwest in summer. Snowpacks were recorded at high elevations between December and May all years: persistently above 700 m a.s.l. and intermittently between 400 and 500 m a.s.l. However, there was a high amount of variation in weather conditions by season, year, location and elevation, which will be described in more detail below.

300
In general, mean annual precipitation increased in a west to east direction, and by elevation, however the rates of change varied due to local conditions. Calvert Island received more precipitation (3353 mm per year) than Hecate Island (2676 mm per year), and for low elevation stations (< 100 m a.s.l.), precipitation increased along a west to east gradient by about 120 mm km -1 on Calvert Island and 150 mm km -1 on Hecate Island (Fig. 3). Precipitation ( Fig. 4). This lapse rate varied seasonally with small differences in temperature between low and high elevation 315 stations during summer months (-0.1 °C per 100 m in August) and large differences in winter (-0.5 °C per 100 m in February). All stations recorded average relative humidities between 80 and 87 %, with high elevation stations reaching lower RH's (e.g. minimum of 16 % daily RH at 'East Buxton') than low elevation stations (e.g. minimum of 37 % daily RH at 'SSN708'). General wind directions were uniform across Calvert and Hecate Islands, except for the areas around Mount Buxton which showed local variability. Highest wind speed was recorded at 'Hecate' station 320 (27.3 m s -1 ), which also recorded highest winds on average (4.4 m s -1 ) and station 'SSN693' was, despite being adjacent to a lake, most sheltered (average of 1.1 m s -1 ). Station specific annual weather parameters are summarized in Table 2.

Temporal variations
Calvert and Hecate Islands had wet (Oct -Apr) and comparatively dry (May -Sep) seasons, with the dry season 335 starting abruptly in May but transitioning gradually into the wet season in September (Fig. 5). Although average monthly rainfall was distinctly lower in the dry season (128 versus 320 mm), rain events still regularly occurred: the longest dry periods never exceeded between 9 and 11 consecutive days, usually in July or August. The wet seasons were marked by large inter-monthly as well as inter-annual variations (58 -524 mm per month), whereas variations were comparatively small during the dry seasons (22 -218 mm per month). November was the wettest, and July the 340 driest month of the year (average of 362 and 78 mm respectively). Frequently recurring high-wind storms (> 10 m s -1 ) from a predominantly southeastern direction characterized the wet season and low northwest as well as southwest winds (< 10 m s -1 ) prevailed during the dry season (Fig. 6). Monthly average air temperature followed an annual cycle ranging between 2.2 °C in February to 14.5 °C in August, with January peaking to 4.3 °C in mid-winter. This seasonal trend was consistent between years, except for the month of February which varied markedly from year to 345 year (between -0.8 and 7.6 °C). Similar to air temperature, relative humidity followed an annual cycle peaking in August-September (90 %) and reaching a low in February (85%), with relatively more year-to-year variability in winter.  (Table 3).

355
which was within the short record period anomalously wet and cool (7.6 °C, 3681 mm, Table 3). However, interpreting these deviations from the climate normal warrants caution, as they are more likely explained by discrepancies in the estimation of climate parameters by ClimateNA. When comparing the observatory's data with projected ClimateNA data for the period between 2016 and 2019 (averaged station specific precipitation and air temperature), they overestimated precipitation at Hecate Island weather stations by 30-700 mm, but underestimated 360 precipitation at Calvert Island by 80-900 mm ( Table 2). The precipitation overestimation on Hecate Island is also apparent when comparing watershed averaged modelled climate normals with measured data, which exceed 1000 mm per year (Table 4)

Streamflow and catchment processes
Average yearly runoff from all watersheds, scaled by watershed area, was 2317 mm, suggesting that ~30 % of 385 precipitation (3246 mm) is not accounted for in surface runoff. The average flux of freshwater from the seven watersheds was 0.1087 km 3 per year. Seasonal runoff patterns followed the wet (Oct -Apr) and dry (May -Sep) patterns of rainfall, with runoff dropping abruptly in May but increasing gradually into the wet season in September reaching highest flows in November and December (Fig. 5). Runoff generated in the wet season accounted for 84 % of total average yearly runoff.

390
Freshwater fluxes were dominated by large storm events: 34 % of all freshwater delivered to the ocean occurred during very high flows that were only exceeded for 5 % of the record (≥ P5 flows), and although none of the rivers fully ceased to flow, baseflows were low with P95 flows approaching zero (Table 4). The majority (92 %) of the ≥ P5 flows occurred during the wet season, with 48 % occurring in November and December. In accordance, 98 % of all days with very low flows (≤ P95) occurred during the dry season of which 52 % occurred during the month of 395 August. Storm events resulted in rapid streamflow responses: average lagtimes (the time from peak rain to the start of rise in streamflow) were less than 8 h for most watersheds; average peaktimes (the start of rise in streamflow to peak discharge) were generally under 12 h (Table 4) respectively. There are potentially three factors that contributed to these differences; first the presence of Mount

415
Buxton with peak elevations of 1012, 680, and 385 m a.s.l. in watershed 703, 693 and 708 respectively, which resulted in elevated precipitation lapse rates on eastern Calvert Island. Second, the majority of large storm events occurred during the wet season when dominant wind direction was from the southeast, and thus there was potential that Hecate Island and western Calvert Island were being affected by a rain shadow from Mount Buxton. Last, the steep relief, exposed bedrock, and sparsely vegetated and shallow soils at high-elevation areas of Mount Buxton 420 contributed to rapid runoff generation and elevated runoff coefficients of 77 and 81 % (versus 60 to 66 % for the other watersheds, Table 4). However, presented watershed averaged precipitation estimates do not include any stations above 740 m a.s.l., and thus are potentially underestimating total precipitation, resulting in overestimated runoff coefficients in watershed 703 and to a lesser extent in 693.
A muted freshet from seasonal snowpacks contributed to stream flow in watersheds 703 and 693. Visual inspection 425 of the hydrograph of watershed 703 showed daily fluctuations of streamflow during dry conditions from April to June, with peak flows occurring ~12 hours after maximum daily air temperature (Fig. 8). These streamflow fluctuations were, compared to peak flows after rainfall events, small (increases of a maximum of 0.2 m 3 s -1 , versus up to 1 to 40 m 3 s -1 during rain events) and could not be distinguished during rainy periods, suggesting snowmelt makes up a small component of annual runoff generation. However, the Mount Buxton snowpack has the potential 430 to contribute to peak flows during rain-on-snow events, especially in watersheds 703, 693 and to a lesser extent in 708.
Based on the general climatic conditions and hydrographs, these watersheds have rainfall-dominated (pluvial) streamflow regimes with peak flows occurring from the early fall through to mid-winter and low flows in the summer dry months (Moore et al., 2012). This contrasts to the freshwater fluxes from one of the major drainage 435 basins within the larger study area, the Wannock river catchment (3900 km 2 , ~50 km east of Calvert Island), which has a snow (nival) and glacial melt dominated regime, with the highest sustained flows in the summer months. Total Wannock river freshwater fluxes reached on average 9.5742 km 3 per year which is around a factor 100 higher than the fluxes from the seven Calvert and Hecate Island watersheds combined (0.1087 km 3 per year). However, yields (scaled by watershed size) were comparable (2317 vs. 2455 mm) (Wannock River daily discharge data, 2021). Despite major differences in flow regimes across the region as described above, Wannock river peak flows generally occurred in the late fall and early winter when atmospheric rivers make landfall (Wannock River annual instantaneous extreme discharge data, 2021), contributing up to 44 % of annual runoff in coastal watersheds (Sharma and Déry, 2019). These atmospheric rivers, combined with warm temperatures and strong winds can melt early season snow and glaciers at higher elevations. It is during such storms that freshwater discharge was 445 synchronized across the region, with corresponding high material flux to the marine ecosystems. It is expected that current yields and streamflow regime will change due to accelerated glacial loss in the Wannock basin (Menounos et al., 2019), a projected increase in the region's mean annual precipitation with less falling as snow (Shanley et al., 2015), and more frequent and intense atmospheric river events (Déry et al., 2009).     Table   2 for station specific PAS model estimates). Rainfall data that were corrected for wind-induced undercatch were on average 12 % greater than the raw data. The largest corrections were applied at 'WSN693_703' (13 %), which was exposed to some of the highest recorded windspeeds (Table 2). In contrast, station 'SSN693' was at a much lower elevation with lower exposure, resulting in overall corrections of 7 %. Some stations exposed to wind required additional anchoring and thus had increased uncertainty due to wind induced tips between 2014 and 2018 ('WSN693_703', 'WSN703_708', 'Hecate'). This was 465 corrected using a tip threshold (> 3 tips per 5 s) to remove faulty tips. However, it was difficult to discern whether some of these tips represented real events even when using adjacent stations as comparison. Following proper anchoring in 2018, 1 to 2 % of the annual rainfall record was flagged as suspect and cleaned. Tipping bucket calibrations indicated there were no shifted values over the course of the study period, but site-specific issues with tipping buckets included blown off lids ('WSN703_708'), blockages from debris ('WSN703' and 'WSN844') and 470 wiring damage by wolves causing false tips ('819_1015'). Resulting issues were addressed and data gaps were filled using linear regression of nearby station records.
The majority of weather station data presented here only measured rainfall (twelve out of fourteen stations, Table 2), and thus there was likely an underestimation of precipitation across the study area. All stations above 400 m a.s.l. measured snow depth to assess the relative contribution of snow, but these data were not used to correct rainfall at 475 these locations. However, according to ClimateNA modelled data, snow made up less than 5 % of annual precipitation at all stations (except 'East Buxton' at 740 m a.s.l.) including those between 400 and 500 m a.s.l. It is therefore likely that any errors associated with not measuring total precipitation at all sites is less than the variation in precipitation across the study area. Snow depth data were corrected by removing obvious spikes, induced through false returns during heavy snowfall 480 and wind, but large data gaps due to sensor failures were not filled with this dataset. It should be noted that the SR50 snow depth sensors had a high rate of failure (1 to 2 years), and thus sensors were replaced annually. In contrast, the air temperature, relative humidity and solar radiation data were without issues. The anemometers occasionally froze in winter causing wind direction records to get 'stuck'. However, these issues were usually short-lived as sensors would melt during daytime.

Uncertainties in the streamflow data
Stream flow measurements are affected by many sources of uncertainty, and thus it is important to both identify and quantify these errors. One potentially major source of error which can be difficult to quantify, is the lack of highflow discharge measurements necessitating the extrapolation of a rating curve from the highest measured flow up to the highest measured stream stage (Coxon et al., 2015;Domeneghetti, 2012). However, the automated salt dilution Channels were stable in watersheds 626, 693 and 1015 for the duration of the data presented, but large storm events resulted in rating curve shifts in watersheds 708, 703, 819 and 844. These watersheds all experienced one to two 495 rating curve shifts within the five-year data record (Fig. 2). However, the automated salt dilution measurement method allowed regular and ample discharge measurements before and after big storms, and thus shifts could be tracked with minimal rating curve data gaps. In addition, peak flows were measured at the upper end of the rating curve in most watersheds. For all watersheds, less than 0.12 % of discharge data points were calculated from the extrapolated part of the rating curve; in other words, less than 2.5 days of the five-year discharge time series are 500 estimated and not measured. On a volumetric basis, between 1 and 3 % of total water discharge from watersheds 626, 693, and 819 was calculated from the extrapolated parts of rating curves, and below 1 % for all other watersheds.
Despite having discharge measurements over the majority of stage measurements in the rating curves, there was significant error in the discharge dataset presented, with mean annual runoff from the seven watersheds potentially 505 being overestimated by 24% and underestimated by 43%. Uncertainty varied by watershed (Table 4), with amounts depending on the spread and measurement uncertainties of individual discharge data points at different positions along the rating curves (affecting the widths of confidence bands). Underlying causes were related to channel stability, turbulence, downstream mixing of salt solution, and general uncertainties introduced by equipment accuracies, which are described in more detail below.

510
Discharge measurement error was always below 20 % at low flows, and generally under 5 %, but occasionally no more than 15 % at high flows using the automated salt-dilution method. The largest factor contributing to the uncertainty of low flow measurements using the velocity-area method, was uncertainty in the velocity readings: despite careful selection of river cross-section sites, ideal sites with minimal flow velocity variation were sparse due to the rivers' steep gradients and complex streambeds. The largest factor contributing to salt-dilution measurement 515 uncertainty was error associated with the salt solution-electrical conductivity calibration factor, which was sensitive to the salt solution's mixing prior to transfer, temperature, and the precision of the calibration equipment (~70 % of total measurement error). The next largest factor was error associated with determining the salt solution dump volume, which due to the system being automated, had to be determined indirectly from pressure measurements (using a pressure transducer inside a stainless-steel collector) and could therefore vary with salinity and within the 520 limits of sensor accuracy. This uncertainty increased with decreasing volume (~30 % of total measurement error).
Further on, discharge measurement error was increased by incomplete mixing of salt solution, or excessive noise in the EC data at downstream sensors. However, most measurements affected were removed from further analysis, except for measurements at the upper end of the rating curves where no other data were available (e.g., high flows at watersheds 703 and 1015 reaching up to 15 % measurement uncertainty).

525
Selecting a stage value to match a discharge measurement for rating curve plotting was difficult when flows were either turbulent and/or the stream level was rapidly rising or falling during the discharge measurement. This error was quantified by propagating the standard deviation of stage values recorded during a discharge measurement to total discharge uncertainty. Fluctuating stage values were mostly a concern for high-flow measurements at https://doi.org/10.5194/essd-2021-423 watershed 819, where stage values could be noisy and rapidly increase or decrease within a 5 to 15 minutes 530 measurement interval, increasing discharge measurement uncertainties by on average 2 % (range of 0 to 9 %). This was also a concern at watershed 844, but to a lesser magnitude (uncertainty increases of 0 to 3 %). All other watersheds were only minorly affected (typically < 1 %).
Rating curve confidence intervals are narrow, resulting in low overall discharge uncertainty, if 1) discharge measurement uncertainties are low, and 2) there is little spread in the rating curve data points. These conditions were 535 met for watersheds 626 and 708, resulting in an overall uncertainty in yearly runoff of -18/-19 % (lower CI) and +23 % (higher CI) ( Table 4, Fig. 2). The rating curves of watersheds 703 and 1015 were of high quality at low to mid flows, but confidence intervals widened above 15 and 1 m 3 s -1 respectively, because of increased measurement uncertainties related to incomplete salt solution mixing at the downstream sites. Measurements taken at watershed 693 all had low uncertainties (< 5 % except at very low flows), but nonetheless a spread in the data (> 3 m 3 s -1 ) 540 increased confidence intervals suggesting there were sources of error unaccounted for in the quantitative uncertainty estimations. Measurements taken at mid to high flows at watersheds 819 and 844 (> 4 and 5.5 m 3 s -1 respectively) had the highest uncertainty due to incomplete salt mixing at the downstream EC sensors and rapidly fluctuating stage during high flows. In addition, unstable channel conditions at watershed 819 caused a rating curve shift once per year, and in some instances the complete rating curve could not be rebuilt before the next shift. All above issues 545 resulted in an overall uncertainty in yearly runoff of -32/-33 % (lower CI) and +74/77 % (higher CI) for watershed 844 and 819 respectively.

Data availability
All discharge and weather data presented in this paper are available at https://doi.org/10.21966/J99C-9C14  and watershed delineations with watershed metrics can be found at http://dx.doi.org/10.21966/1.15311 550 (Gonzales and Giesbrecht, 2015). Both data packages include a readme file detailing data file content and contact information for further details.

Code availability
Calculations were done in R: the rating curve script is available at https://github.com/HakaiInstitute/RatingCurve and the scripts used for weather data QC can be found at https://github.com/HakaiInstitute/wx-tools.

8 Conclusions
The hydrometeorological data of the Kwakshua Watersheds Observatory fills a critical data gap on the outer coast of the northeast Pacific coastal temperate rainforest (NPCTR) of North America and over time will help to better understand the hydrology of the region. The hydrometeorological dataset presented is part of a greater interdisciplinary effort to better understand this region, including oceanography, biodiversity, ecology, climate and 560 watershed processes. High quality data are assured using state-of-the-art technology and thorough data quality procedures. Discharge data presented here are unique because the novel method of measuring streamflow accounts for measurement error from various sources, which is often not quantified in publicly available datasets. The dataset https://doi.org/10.5194/essd-2021-423

590
The authors declare no competing interests.

Acknowledgements
We gratefully acknowledge that this work was conducted on the traditional, ancestral and unceded territories of the