Initial Distribution and Interannual Decrease of Suspended Sediment in a Two‐Basin Lake Following a Massive Mine Tailings Spill: Quesnel Lake, BC, Canada

On 4 August 2014, a tailings dam failure at Mount Polley Mine in central British Columbia led to the largest short‐term release of mine waste into a lake ever recorded. Once released by the breach, slurry entered the smaller of Quesnel Lake’s two basins, called the West Basin, from which the lake drains northwest into the Quesnel River. An estimated 38,000 ± 11,000 tonnes of fine solids remained suspended in the West Basin on 10 September 2014, 37 days postspill. This decreased to within background levels (<300 tonnes) by June 2015, by which time ∼4,000 tonnes of sediment had flowed from the West Basin into the Quesnel River, and ∼31,000 tonnes had been transported east into the main basin, the direction opposed to the mean hydraulic gradient. Here, we evaluate sediment transport in Quesnel Lake following the Mount Polley tailings dam spill. We apply conservation of mass in two ways: using data collected between 12 August 2014 and 15 October 2020 to estimate suspended sediment mass and mass flows into and out of the West Basin; and using one‐ and two‐basin completely mixed models. Suspended sediment concentrations were highly elevated through the first three seasons postspill, then fluctuated within a gradually decreasing, seasonal cycle. The observed mass trend and an analytical mass balance model of a simplified, two‐basin system suggest that suspended sediment mass will fall to below the median prespill mass by 7 ± 5 years postspill.

and never found), and releasing 12 million m 3 of toxic, metal-laden tailings which deposited in the Paraopeba River (Vergilio et al., 2020). By released volume, the first of these two disasters is the largest recorded to date; the latter is the third largest yet (WISE, 2021). Rivers were the receiving waters affected in both cases (do Carmo et al., 2017;Vergilio et al., 2020).
The second largest mine tailings spill on record affected two lakes. In the early hours of 4 August 2014, a breach formed in the perimeter embankment dam surrounding the tailings storage facility at Mount Polley Mine, a copper and gold producer in central British Columbia, Canada. Fortunately, no mine personnel were on or below the affected section of the dam when it failed, and the event did not result in loss of human life. Over the course of approximately 12 hr, however, the breach released 7.3 million m 3 of tailings solids, 10.6 million m 3 of supernatant (water-cap), and 6.5 million m 3 of interstitial water into the surrounding environment (MPMC, 2015). Of the ∼25 million m 3 that escaped through the breached dam, about one quarter (6.4 ± 0.5 million m 3 ) went into nearby Polley Lake and was held behind a broad, geotechnically unstable pile of tailings called the "Polley Plug" which blocked the lake's natural outflow into Hazeltine Creek (MPMC, 2015). The remainder of the flood reshaped the ∼9-km-long Hazeltine Creek channel and discharged ∼21 million m 3 into Quesnel Lake (Petticrew et al., 2015). Offshore of the Hazeltine Delta, the downward-flowing debris eroded two subaqueous channels into the sloped lake bottom (Kostaschuk et al., 2021). Solids from these debris flows accumulated in a layer up to 15 m thick over deeper, flat-bottomed regions in the West Basin, Quesnel Lake's smaller, downstream basin (Kostaschuk et al., 2021). This bottom deposit, which accounted for ∼99.7% of the total mass of solids that the spill introduced to Quesnel Lake, was predominantly composed of silt-sized (2-62.5 μm) particles (Tetra Tech EBA, 2015).
The main interest of our study was the fate of the remaining 0.3%-roughly 30,000 tonnes of predominantly claysized (<2 μm) particles-which within days of the spill were suspended in a plume that had spread throughout the hypolimnion below 30 m depth in the West Basin (Petticrew et al., 2015). (Note that from here on, sediment mass is given in units of megagrams; 1 Mg = 1,000 kg = 1 tonne). Quesnel Lake was strongly stratified, as it always is in early August, with its surface temperature near 20°C through a ∼10-m-deep epilimnion, its thermocline at ∼15 m below the surface, and its hypolimnion near 4°C below 25 m depth. In the West Basin's hypolimnetic plume, the physical properties of spill-contaminated water were distinct from those of deep water in the rest of Quesnel Lake in several ways: the plume was warmer, at 6°C-7.5°C; its turbidity ranged from 200 to >1,000 NTU, compared to a natural background level of 0.4 NTU; and its conductivity was 160 μS/cm, well above the background of 110 μS/cm (Petticrew et al., 2015). Suspended sediment in water samples collected from 30, 60, and 90 m depth in the West Basin on 18 September 2021, 45 days postspill, had a median grain size (D 50 ) slightly smaller than 1 μm, a D 10 around 0.5 μm, and a D 90 between 3 and 7 μm (Petticrew et al., 2015).
A 35-m-deep sill separates the West Basin from the main basin, providing a 10 m vertical gap through which spill-contaminated hypolimnetic water could pass into the main basin. In the months following the spill, a bifurcated plume propagated eastward at ∼1 cm/s, ∼20 m thick below the thermocline, and ∼30 m thick at the lake bottom (Petticrew et al., 2015). Strong currents are known to exchange water over the sill between the West Basin and the main basin during periods of sustained wind forcing (Laval et al., 2008); these bidirectional exchange flows would have contributed to suspended sediment transport over the sill.
Considering the scale of the Mount Polley tailings spill and the notable impact it wrought upon the West Basin through late 2014 and into early 2015, it is remarkable that Quesnel Lake's temperature, turbidity, and conductivity had returned to near prespill conditions by summer 2015 (Hamilton et al., 2020). That said, turbidity in the West Basin hypolimnion was anomalously high (2-2.5 formazin turbidity units [FTU]) each autumn and spring of 2015-2017 compared to prespill levels. And, in summer 2016, a several centimeter-thick, cloudy, mobile layer was observed overlaying more consolidated material in cores taken from the heavily impacted bottom region of the West Basin, but not cores taken from the bottom of the main basin (Hamilton et al., 2020). Commonly referred to as "benthic boundary layers" or "nepheloid layers," and typically found overlying bottom sediments in the deepest parts of a lake, such layers tend to have elevated turbidity due to particles kept in suspension by turbulence (Gloor et al., 1994). Hamilton et al. (2020) showed that increased autumn turbidity was associated with baroclinic seiching and hypothesized that the source of this turbidity was the cloudy, mobile layer observed in the West Basin.
This study is meant to extend earlier observations of turbidity made by Petticrew et al. (2015) and Hamilton et al. (2020) following the Mount Polley tailings dam spill. But whereas those earlier studies also considered changes in temperature and conductivity associated with the 2014 spill and subsequent effluent discharge, this study considered only turbidity-or more specifically, the mass concentration of suspended sediment estimated from turbidity using a site-specific relation. While changes in conductivity are noteworthy, as they imply higher dissolved solid loading to Quesnel Lake, chemical changes are not the focus of this study. Our work is meant to complement earlier publications on Quesnel Lake with a deeper, more quantitative analysis of suspended sediment. To that end, we conducted a sediment mass balance based on field data observations from the first 9 months postspill; quantified (in terms of sediment mass) the seasonal variations in turbidity which Hamilton et al. (2020) observed beyond 10 months postspill; established what the background level of naturally occurring suspended sediment is in Quesnel Lake; and assessed suspended sediment assimilation for the first decade postspill using mathematical models in order to predict when suspended sediment concentration will return to background levels.
Quesnel Lake is currently used as the receiving water body for effluent discharge by Mount Polley Mining Corporation (MPMC), as regulated under the Province of British Columbia's Environmental Management Act (permit no. 11678). For the Quesnel River watershed, ongoing mine water management decisions would benefit from a postspill assessment of suspended sediment mass, years 1-6 (2014-2020). More broadly, the science of lake sediment transport is missing a quantitative description of what happens when a large amount of sediment goes into a two-basin, deep fjord lake.

Study Site
Quesnel Lake (volume: V = 41.8 km 3 ; surface area: A = 266 km 2 ), shown in Figure 1, extends from the foothills of the Cariboo Plateau into the Cariboo Mountains of central British Columbia, Canada. Viewed on a map oriented east-up, the lake's three arms resemble a lowercase "y." Its deepest point is in the East Arm, 512 m beneath the surface of the lake (surface elevation: 725 m above sea level), making it the deepest fjord lake in the world and the third deepest lake in North America (Gilbert & Desloges, 2012). Like most fjord lakes, it is long (east-west thalweg length: 95 km), narrow (mean width: 2.7 km), and has multiple basins (two). The West Basin (V d = 1 km 3 , A d = 22 km 2 ) is the hook at the bottom of the lowercase "y" (Figure 1, main panel). It curves from a westerly to a northerly orientation over 20 km, broadening and deepening from a 35 m sill near Cariboo Island to its deepest section offshore of Hazeltine Creek (Figure 1), then becoming shallow and narrow toward the outflow. Its maximum depth prior to the spill was 121 m (Gilbert & Desloges, 2012); the 4 August 2014 spill entered the West Basin via Hazeltine Creek, and reduced the basin's maximum depth to 108 m (MPMC, 2015). Near the town of Likely in the northernmost reach of the West Basin, Quesnel Lake drains into the Quesnel River (study period mean annual flow = 140 m 3 s −1 ; ECCC, 2020a, 2020b).
Like most temperate lakes, Quesnel Lake is classified as dimictic, although the deepest parts of the lake may be better described as cold monomictic (Laval et al., 2012). The West Basin mixes through its full depth each autumn and spring; taken together, these two periods of mixis make up around 20% of a typical yearly cycle (Granger, 2020). The heat content of the lake reaches a maximum in late summer, at which time a temperature difference of as much as 10°C can be measured between 10 and 20 m depth (Laval et al., 2012). The autumn breakdown of stratification coincides with strong wind forcing, resulting in large amplitude baroclinic motions (Thompson et al., 2020). As far as we are aware, there are no records of the East Arm ever having frozen over completely. Quesnel Lake does receive partial ice cover during some winters, with winter ice cover occurring more commonly in the West Basin than elsewhere in the lake (Hamilton et al., 2020).

Field Data
In the wake of the August 2014 spill, several organizations either initiated or increased water quality monitoring in Quesnel Lake and the Quesnel River. This study used data collected by provincial and federal governments, industry, and public research institutions (Table 1); the authors are affiliated with a collaboration that we herein refer to as the Environmental Damages Fund (EDF), named for the organization which provided our main financial support through Environment and Climate Change Canada (ECCC).
EDF data include vertical profiles collected in Quesnel Lake using three Seabird Electronics conductivity, temperature, depth (CTD) profilers, all equipped with Seapoint turbidity sensors (Table 1). Our West Basin prespill CTD data set consisted of 16 profiles from 12 dates between 2006 and 2012, with data from spring through early autumn, and no data from late autumn or winter (note: throughout the text, the terms summer, autumn, winter, and spring correspond to the solar year, except when followed by a state of stratification, e.g., spring mixis or winter inverse stratification). The 10 CTD stations shown in Figure 1 were profiled on a biweekly basis during the late April to early December field seasons; here we include 100 of these transects from 10 September 2014 onwards, the date of the first available full-depth West Basin profiles. The prespill data and the early years (2014)(2015)(2016)(2017) of postspill data were collected using an SBE19plus (SN 4057); the 2018 CTD data using an SBE19plusV2 (SN 7035); and the 2019-2020 CTD data using an SBE19plusV2 (SN 7098). Turbidity (Tu, FTU) data from each CTD have been corrected to formazin standard. Lake with 20 m depth contours. Inset (i) shows the location of the lake in British Columbia, inset (ii) shows the location of the West Basin in the lake, and inset (iii) is an along-thalweg depth plot for the West Basin indicating locations of conductivity, temperature, depth (CTD) stations and moorings (also indicated in main panel). Note that ST11 and ST7 have shallow cast depths indicated by dotted lines in (iii), as these two are located near shore.
Throughout the first winter postspill, MPMC used a CTD profiler (YSI EXO2 Sonde) to measure physical water properties in Quesnel Lake. As far as we are aware, these profiles were the only high vertical resolution turbidity data collected during this period. Herein, we include turbidity data from 12 profiles collected between 1 November 2014 and 30 May 2015 at station QUL-66a, and 1 profile collected on 10 December 2014 at station QUL-112a ( Figure 1). Station QUL-66a turbidity profiles compare well to EDF turbidity profiles (SBE 19plus,SN 4057; formazin-standardized values) taken within a few days from EDF stations ST9 and ST10 in November 2014 and May 2015.
Continuous, in-lake turbidity data came from EDF moored sensors (two JFE Avanatech Co. ACLW-USB, and two RBR virtuoso with Zebra-Tech wipers; Table 1). Between November 2014 and October 2019, the JFE turbidity sensors were deployed in the West Basin (4 and 38 m depth at mooring M3, Figure 1). In October 2019, the JFE sensors were moved to 28 and 53 m depth at mooring MTu. In November 2019, the RBR sensors were added to mooring MTu at depths of 3 and 78 m. We needed to standardize moored turbidity data to CTD turbidity data, as each moored sensor exhibited a unique response to changes in turbidity. To do so, we used in situ CTD turbidity data from the biweekly transects. First, we identified data points for comparison: CTD data from the station and depth nearest each moored sensor, along with the corresponding moored sensor data from the date and time of the CTD profile. We then calculated a slope and offset for each moored sensor using linear, least squares regression; and finally applied these corrections to the moored sensors' turbidity time series.
ECCC monitors flow (Q, m 3 s −1 ) in the Quesnel River through the Water Survey of Canada (WSC), and turbidity through the Federal-Provincial Freshwater Quality Monitoring and Surveillance Network (FWQMS). Historical, daily averaged Q was available for the Quesnel River at Likely Bridge (WSC station number 08KH001) for the years 1925-2019 (ECCC, 2020a); with provisional, real-time data available for 2020 (ECCC, 2020b). FWQMS samples were collected from the surface of the Quesnel River at Likely Bridge weekly in the year following the spill, and then biweekly or monthly thereafter, with turbidity data available to March 2020. At the Quesnel River Research Centre (QRRC), 1.5 km downstream of Likely Bridge (Figure 1), MPMC installed a YSI EXO Sonde on 12 August 2014 which recorded turbidity every 15 minutes (MPMC, 2015). EDF began water sample collection at the QRRC on 26 August 2014, and measured inorganic suspended particulate matter (SPM, g m −3 ) by filtering samples and drying filters at a temperature between 60°C and 90°C. SPM data were available to 15 October 2020.
Bathymetric data were available from a sonar survey made prior to the August 2014 spill, and consisted of 3118 depth points in Quesnel Lake, of which 550 were in the West Basin. To obtain a hypsometric curve for the West Basin, which we needed for our sediment mass balance, we generated a grid with 10 m horizontal and 1 m vertical resolution by interpolating between these 550 points. We then adjusted the hypsometric curve to reflect how the spill reduced the depth of the West Basin. To match the vertical resolution of our hypsometric curve, we averaged CTD turbidity data into 1 m depth bins.

Sediment Mass Balance
Being oligotrophic, Quesnel Lake belongs to a class of lakes which contain predominantly minerogenic particles, such as silts and clays (Håkanson & Jansson, 1983); with less abundant organic (i.e., low loss on ignition) and inorganic, biogenic particles (Wetzel, 2001). That said, minerogenic matter is not the only source of turbidity in Quesnel Lake. A vertical profile collected on 3 June 2008 near station ST9 in the West Basin ( Figure 1) showed elevated turbidity (plotted in Figure 3a as c, g m −3 ) between 1 and 3 m depth, with chlorophyll concentration likewise elevated at the surface (∼1 μg L −1 above a background of ∼0.25 μg L −1 , not shown), which indicates that higher surface turbidity was partly due to algal growth. For our mass balance, we have assumed that suspended organic matter had a negligible contribution to turbidity for two reasons: (a) suspended sediment from the spill was present at a much higher c than the highest observed for chlorophyll and (b) the 9-month period of interest falls mostly within those seasons (autumn, winter, and spring) that tend to have less abundant biological activity.
Earlier studies examined the changes in Quesnel Lake's turbidity that followed the Mount Polley tailings dam spill to establish that suspended, fine sediment spread through the hypolimnion of the West Basin within days, entered the main basin within weeks, and reached the Quesnel River within months (Hamilton et al., 2020;Petticrew et al., 2015). Here, we have gone a step further by using turbidity, some of which appeared in these earlier studies, to quantify sediment mass transport from the West Basin to the main basin and the Quesnel River. Making the step from qualitative to quantitative observations meant converting turbidity, an optical property of water, to suspended sediment mass concentration, a physical one. This conversion needed to be site specific, as the relationship of turbidity to mass concentration varies for differing sediment compositions or particle size distributions (Jastram et al., 2010).
Estimating c from Tu is a well-studied topic due to its wide range of applications, including calculating stream loads, which we have done here for the Quesnel River. General interest has led to the development of both univariate (based only on Tu) and multivariate (based on Tu and other parameters) approaches (Jastram et al., 2010).
Here, we employed a univariate approach, estimating c from Tu as with offset taken as the turbidity of distilled water (Tu ∅ = 0.1 FTU), and slope (k c , g m −3 FTU −1 ) from linear, least squares regressions of Tu to SPM data. Note that while both SPM and c have units of g m −3 , SPM is defined more specifically as the concentration of matter which is retained on a 0.45-μm pore size filter paper. For our site-specific c(Tu) relation, we used the lower bound, expected, and upper bound k c (0.66, 0.78, and 0.93 g m −3 FTU −1 , respectively) that Granger (2020) developed in the laboratory using sediment which had been collected in July 2016, with a slow-corer, from the deep part of the West Basin (a plot of these laboratory data may be found on p. 20 of that document). Tetra Tech EBA (2015) similarly gave a lower bound estimate of k c = 0.66 g m −3 FTU −1 for this sediment based on in situ water column stability.
We defined suspended sediment mass (m) locally for the water column of the West Basin, calculating upper and lower bounds of m from vertical turbidity profiles collected at five CTD stations (ST6, ST7, ST9, ST10, and ST11, Figure 1) as follows: from each set of five profiles, make an upper bound vertical Tu profile composed of the maximum observed values across all profiles for each 1 m depth bin, and likewise use minimum Tu to make a lower bound vertical profile; add/subtract sensor error to upper/lower bound Tu profiles; convert these to vertical c profiles using Equation 1 with upper/lower k c estimates; and integrate each through the West Basin's area-depth curve (A(z), z axis origin is at the surface and is positive downwards, z b = 108 m): where c n and V n are the sediment concentration and volume of the nth depth bin, respectively. Our best estimate (expected m value) was based on the median turbidity across all profiles for each bin, with c(z) from Equation 1, the same A(z) curve as we use for upper and lower m bounds, and k c = 0.78 g m −3 FTU −1 . In the few cases where we have calculated m by Equation 2 for natural sediment, we have stated k c explicitly.
Conservation of mass gives that the rate of change of sediment in a basin is equal to the sum of mass flows across its boundaries (due to the physical processes of loading, settling, and flushing) and chemical transformation within the basin. In Section 3, we evaluated mass flows associated with three boundaries of the West Basin during the first 9 months postspill: to the north, into the Quesnel River ( , Mg day −1 ); along the bottom, due to settling (̇) ; and to the east, into the main basin of Quesnel Lake (̇) . For suspended sediment mass flow in the Quesnel River over these 9 months, we have used historical, daily averaged flow (ECCC, 2020a) and MPMC turbidity data collected at the QRRC to estimate sediment mass flow with c from Equation 1 using k c = 0.78 g m −3 FTU −1 . We have assumed that MPMC turbidity represents an averaged c across the width and through the depth of the channel, as the river is well mixed by several sets of rapids between the outflow and the QRRC.
Suspended sediment collected in bottle samples from 30, 60, and 90 m depth in the West Basin on 18 September 2014 had a median particle diameter slightly smaller than 1 μm, with a density slightly greater than ρ p = 2,500 kg m −3 (Petticrew et al., 2015). For such a particle, the Stokes' settling velocity (w s , m s −1 ) is where g = 9.8 m s −2 is the gravitational constant, ρ f = 1,000 kg m −3 is fluid density, and μ = 0.0015 kg m −1 s −1 is the dynamic viscosity of water (Håkanson & Jansson, 1983). Using this velocity, which assumes spherical particles, no interaction between particles, and a quiescent environment, we calculated an approximate settling mass flow rate (̇) for the date of each West Basin c(z) profile aṡ where L(z) is the isobath length at depth z, and A n is the bottom area of the nth depth bin. For dates between c(z) profiles, we estimated using logarithmic interpolation; that is, linearly interpolating in ln (̇) . We do not give upper and lower bounds for ; rather, we treat these as order-of-magnitude estimates, reflecting their uncertainty. The exact conditions affecting settling are complex, as flocculation increases settling velocity, and currents may enhance or negate a particle's downward motion through the water (Gloor et al., 1994;Mehta, 2014;Scheu et al., 2015).
We estimated suspended sediment mass (Mg) transport out of the West Basin water column for the initial months following the spill from the mass flows , , and (Mg day −1 ). In the first and second case, we assume a constant rate of mass flow for each day, multiply from Equation 3 or from Equation 5 by Δt = 1 day, and sum over the number of days in each period. In the third case, we apply the mass balance: to assess cumulative mass flow eastward from the West Basin over the Cariboo Island sill into the main basin of Quesnel Lake by residual. We take the change in observed sediment m (Table 2) for each period, and subtract cumulative and to get cumulative (Table 3).

Mathematical Models
The severity and persistence of environmental impacts resulting from a sudden input of sediment depend largely on the assimilative capacity of the receiving waters-their ability to flush away, bury, or absorb sediment (Chapra, 1997). We have used continuously stirred tank reactor (CSTR) models to investigate the interannual trend in m. Such models use constant terms, and so cannot predict seasonal variations. However, because the lake being modeled mixes from top to bottom at least once per year, we assume that these models can predict the net effect of seasonal processes.
We noted in Section 2.1 that Quesnel Lake is seasonally stratified-incompletely mixed through most of its yearly cycle-and so any completely mixed model is an imperfect representation of the system. Here, we assume that for time scales of several years or longer, provided the lake mixes completely at least once per year, a CSTR model can capture the net effect of seasonal processes and still inform our understanding of the long-term trend in suspended sediment mass.
The sediment mass balance for a single, completely mixed basin with constant inflow equal to outflow (Figure 2a) is given by  Note. Error estimates are given for mass based on EDF CTD transect data, as these take into account spatial variation in suspended sediment concentration based on five casts and include upper and lower bound c(Tu) correlations (Section 2.3). The remaining mass estimates are based on k c = 0.78 g m −3 FTU −1 unless otherwise indicated. The description "summer minimum" refers to the yearly minimum which occurred during the summer stratified period, some years in the solar summer, others in autumn. The bold entries indicate the event which the article is about. a k c = 1.3 g m −3 FTU −1 . b k c = 1.2 g m −3 FTU −1 . where W(t) is loading (g s −1 ), and r is a first-order reaction constant (s −1 ) for chemical transformation of solid m into soluble m. Grouping constant terms, we have with the eigenvalue λ, or characteristic value of mass decay for sediment in the system, given by where V/A = H, average depth (m). The inverse of λ is the e-folding time, for 63% removal (1 − e −1 ≈ 0.63) of mass. For dissolved matter, Equation 9 retains the flushing term (Q/V) but not the settling term (w s /H), and provided the substance is not removed by chemical reactions (r = 0), its e-folding time will be the basin's residence time (t r = V/Q). We consider only settling and flushing in our analysis, as Quesnel Lake water has a pH of 7.4 ± 0.2 (Laval et al., 2012), and the sediment introduced by the Mount Polley tailings dam spill is considered to be chemically stable in this range (SRK Consulting Inc., 2015).
In the case of sudden loading, W(t) can be mathematically represented in terms of the Dirac delta function (δ(t), s −1 ): where m i is the mass of sediment released into the basin. The solution to Equation 8 becomes For a feedback system of two basins (Figure 2b), flow between the basins is bidirectional, with exchange flow (Q x ) from the downstream (denoted by subscript d) to the upstream (subscript u) basin, and exchange flow plus river flow (Q x + Q) in the opposite direction. This yields a system of two linear, homogeneous, ordinary differential equations: with coefficients (α) given by the transport terms:   The general solution to this two-basin system has two decay terms, denoted by the subscripts f and s, which correspond to the fast and slow system response: with eigenvalues λ f and λ s given by and eigenvectors η f and η s given by Having calculated λ f,s and η f,s , one specifies the initial concentrations c u (0) and c d (0) and solves Equation 19 for k f and k s .

Observations
On 10  West Basin surface c gradually increased through autumn 2014 until, by November, the initially clear surface water had become visibly green (Figure 4a; see Hamilton et al. [2020] for satellite imagery and an account of surface mixed layer deepening and entrainment of the turbid, hypolimnetic water that led to this "greening"). At 93 days postspill, suspended sediment m in the West Basin's water column was about a quarter of what it had been two months prior (9,200 ± 1,900 Mg on 5 November 2014; Table 2). By 10 December, with 7,600 Mg remaining in the West Basin's water column, c was near 7 g m −3 through depth near CTD station ST9 (Figure 3b), consistent with sediment c near 7 g m −3 in the Quesnel River (FWQMS data, Figure 4c). A steady decrease in West Basin sediment c through the first winter postspill was interrupted by a brief period of increase at the beginning of spring 2015 (Figure 4a). By the end of spring mixis in late April 2015, sediment m had fallen to ∼1,000 Mg (Table 2). Sediment m continued to decrease through summer 2015, eventually reaching a late-August minimum of 140 ± 30 Mg (Figure 4b and Table 2). Thus, in the year following the spill, West Basin suspended sediment m underwent a 2 order of magnitude decrease (Table 2).
Suspended sediment mass flow out of the West Basin and into the Quesnel River between 10 September 2014 and 20 March 2015 accounted for 4,000 Mg, or just over 10% of the total decrease in West Basin m (Table 3). A smaller fraction of the initial West Basin m decrease can be attributed to settling (Table 3) (Figure 3e). A physical mechanism for transporting sediment over the Cariboo Island sill can be seen in acoustic Doppler current profiler data from that location (Figure 8a in Hamilton et al. [2020]). Oscillating, bidirectional exchange flow (seiche pumping), first documented by Laval et al. (2012), was observed through late 2014 and into spring 2015. Turbulent mixing between the hypolimnion and epilimnion is likely to occur where flow is most constricted, at the sill, and the net effect of the back-and-forth motion of seiche pumping is sediment transport from areas of higher c to areas of lower c: a strong sediment c gradient, from high c in the West Basin to low c in the main basin, persisted until spring 2015, consistent with net eastward for the period (Table 3).
Through the rest of the study period, summer 2015 to autumn 2020, we saw seasonal fluctuations in sediment c and m which followed the annual progression of stratification and mixis in the West Basin (see Hamilton et al. [2020] for detailed observations of turbidity up to autumn 2017). In 2016 and 2017, we recorded annual maximum sediment m in autumn; in 2018, 2019, and 2020, maximum m occurred in spring (Table 2). We did not have vertical turbidity profiles (CTD data) from periods of winter inverse stratification, but data from moored turbidity sensors indicate that sediment c decreases slightly between autumn and spring mixis. Sediment m would have reached a local minimum at some point each winter, though not as low as the summer minimum due to the shorter duration and weaker stratification of winter. Thompson et al. (2020) note that the wind over Quesnel Lake is at its most calm during August, and particle settling through a more quiescent water column may have contributed to annual m minima each summer; conversely, wind forcing was found to be strongest in April and November, around those times of year when we observed sediment m maxima (Table 2).
Nearing the end of the study period, in spring 2020, point samples of Quesnel River inorganic SPM and continuous West Basin turbidity data provided us an estimate of the linear c(Tu) coefficient, k c , for natural sediment, as well as upper bound estimates for natural sediment c in Quesnel Lake and sediment m in the West Basin. Note that the main basin receives sediment from three large rivers (the Horsefly, Mitchell, and Niagara; Figure 1). The hydraulic gradient causes sediment to be drawn into the West Basin, which we observed in July 2020. Between 2015 and 2019, Quesnel River flow peaked in either late May or early June during spring freshet. A deep snowpack brought a strong freshet in spring 2020 which was followed by heavy precipitation in early summer. This led to record high water levels and increased terrestrial sediment loading to Quesnel Lake, and a 1-in-30 years flow in the Quesnel River (Figures 5a-5c). Peak river flow (Q = 585 m 3 s −1 ) occurred on 8 July 2020, coinciding closely with maximum sediment c in the West Basin epilimnion and the Quesnel River (c = 5.5 g m −3 ). CTD data collected on 17 July 2020 show higher turbidity at the Cariboo Island sill (station ST5) than in the West Basin, indicating westward due to sediment transport by hydraulic flow within the epilimnion (Figure 5a). We expect the main source of this sediment to have been the Horsefly River, as sediment plumes originating from the more distant Mitchell River and Niagara Creek have not been recorded as far west as the Cariboo Island sill (Hamilton et al., 2020).
Natural sediment would have been the main contributor to turbidity in July 2020, 6 years after the spill. At peak river flow on 8 July, the hydraulic residence time of the West Basin epilimnion was on the order of 10 days, which suggests that most of this sediment went into the Quesnel River, while relatively little stayed behind in the West Basin. The blue dashed line in Figure 5d shows the correlation of West Basin surface turbidity to Quesnel River inorganic SPM data for the period 21 November 2019 to 14 September 2020, k c = 1.3 g m −3 FTU −1roughly double that of the lower bound k c value for spilled sediment, and well above that of the upper bound (0.93 g m −3 FTU −1 ). Note that if we instead used total SPM (inorganic plus organic) for the same period, this correlation gives k c = 1.4 g m −3 FTU −1 . Using k c = 1.3 g m −3 FTU −1 in Equation 1, we apply Equation 2 to the horizontally averaged West Basin c(z) profile from 17 July 2020 to get m = 900 Mg, our upper bound estimate (from a 1-in-30 years event) for the naturally occurring range of sediment m (Table 2).
CTD turbidity data collected prior to the 4 August 2014 spill provide our median West Basin m estimate, used in Section 3.2 to assess recovery time. The West Basin typically showed low turbidity (<0.4 FTU), with some instances of higher turbidity (1.5 FTU) near the surface during spring (Hamilton et al., 2020). Based on k c = 1.3 g m −3 FTU −1 , hypolimnetic sediment c was typically <0.2 g m −3 in the West Basin (Figure 3a). In the main basin, hypolimnetic sediment c was similarly low, with riverine sediment input bringing slightly higher c to the deeper water (Figure 3d). Using the median c observed in the prespill data at each depth in the West Basin, Equation 2 gives an approximate prespill median m of 130 Mg (Table 2), which divided by V d gives c = 0.13 g m −3 .

Modeling
The suspended sediment mass time series (Figure 4b) shows that through the seasonal phase (mid-2015 onward), suspended sediment concentrations in the West Basin of Quesnel Lake gradually returned toward their prespill median. In the two-basin model, we described in Section 2.4, terms for fast and slow decay are given by Equation 19; we expect the latter term to be important for interannual decrease. Before we discuss the parameters needed for a two-basin model of Quesnel Lake, those needed to obtain a particular solution to Equation 19, we will begin with the simpler, one-basin model for comparison (Figure 2a). Consider the solution for impulse loading to a CSTR, given by Equation 12, in which flushing (Q/V) and settling (w s /H) together make up the eigenvalue λ. Flushing is the inverse of Quesnel Lake's residence time ( −1 = 9.4 −1 a −1 ). A settling velocity of w s = 5 cm day −1 (Section 2.3) divided by a mean depth of 157 m gives a settling term of similar scale to flushing (w s /H = 8.6 −1 a −1 ), resulting in an eigenvalue λ = −4.5 −1 a −1 . This means that Quesnel Lake's e-folding time or "spilled sediment residence time" is around 5 years (i.e., it would take 5 years to decrease from a mass of m i to a mass of e −1 m i ≈ 0.37 m i ).
In the two-basin model of Quesnel Lake, the main basin is upstream, while the West Basin is downstream. The upstream basin has a significantly greater volume (V u = 40.8 km 3 , V d = 1 km 3 ) and surface area (A u = 243 km 2 , A d = 23 km 2 ). River flow is unchanged from the previous case, equal to the mean annual flow of the Quesnel River (Q = 140 m 3 s −1 ). Settling velocity will also remain unchanged at w s = 5 cm day −1 . To parameterize exchange flow, we can use the result of Laval et al. (2008) that during summer stratification, the West Basin hypolimnetic residence time is 6-8 weeks due to seiche driven exchange flow. The volume of the West Basin hypolimnion in late summer is roughly half that of the basin, and a 6-8-week hypolimnetic residence time equates to an exchange flow of comparable magnitude to river flow (Q x = Q).
We initialize the two-basin model for t (0)  The particular solution for c d (t) is plotted in Figures 4b and 4c. As with the two phases observed in the West Basin following the spill, the solution predicts an initial period of rapidly decreasing c transitioning into a slow, sustained decrease. Because the model uses constant transport terms, it cannot capture the seasonality of the real system. Nonetheless, three key aspects of the real system are reflected in the model: the decay rate of the initial phase (given by the eigenvalue λ f ); the timing of the transition between phases; and the decay rate of the seasonal phase (given by λ s ). The similarity between the two-basin model's slow decay and the observed interannual trend suggests that suspended sediment from the spill is being assimilated (flushed away or removed by settling) by both basins of Quesnel Lake, rather than its effects having been confined to the West Basin.
This two-basin analytical model uses seven parameters: Thus, for a spill of such fine sediment into Quesnel Lake, exchange flow mainly affects the short-term system response, while particle settling mainly affects the long term. There is a straightforward physical explanation for why varying Q x has a negligible effect on λ s . Exchange flow is not a sink; it relocates but does not remove sediment from the lake. The transition from fast to slow decay occurs once sediment c has been equalized between the upstream and downstream basins by exchange flow. For the base case, as well as for Q x = Q ± 116 m 3 s −1 , the slow eigenvalue is nearly identical to the eigenvalue in the one-basin model (λ = −4.5 −1 a −1 ). The third and fourth cases imply that the parameterization of w s is important when evaluating the long-term trend. We used a Stokes' settling velocity of 5 cm day −1 for the modeled sediment m shown in Figure 4b and found a good agreement between the slow decay predicted by the two-basin model and the trend shown in the CTD data.
As noted, suspended sediment c was typically very low in the West Basin before the spill (median c = 0.13 g m −3 ). Once spilled sediment c falls below this median prespill c, we consider the lake to have returned to background suspended sediment levels. The one-and two-basin models we have described give similar long-term decay rates, so either may be used to assess how long this return will take. In the case of the one-basin model, we solved Equation 12 for the time at which c(t) = 0.13 g m −3 . Taking m i = 38,000 Mg (the sediment mass estimated to be in the West Basin on 10 September 2014), and V = 41.8 km 3 , the one-basin model predicts that spilled sediment c will fall below the median prespill c after one decade postspill (actual time using these parameters: 9 years, or twice the e-folding time). The two-basin model predicts a somewhat shorter return to background c of just over 7 years, using a settling velocity of 5 cm day −1 . If we instead use w s = 42 cm day −1 , the two-basin model gives a much shorter time of 1 year and 3 months, while w s = 1.2 cm day −1 gives 12 years. Based on these results, we have concluded that our turbidity-based method will no longer be able to detect the fine sediment released into Quesnel Lake during the 2014 Mount Polley tailings dam spill by 7 ± 5 years after the event.

Discussion and Concluding Remarks
In Section 1, we mentioned that a nepheloid layer had been observed in summer 2016, in the deep part of the West Basin, overlaying bottom sediments where the tailings deposit from the spill is several meters thick. Particles in this nepheloid layer remained in suspension near the bottom, so for our two-basin model they would be counted as part of the downstream basin's sediment mass. But this mass was contained in a 10-15-cm-thick layer, which we would not have detected in vertical turbidity profiles collected with our CTDs, because the turbidity sensor on each CTDs is ∼40 cm above the bottom of the instrument. The observed suspended sediment m falls well below modeled m each summer of 2015-2020 (Figure 4b), partly because we were unable to detect sediment m in the nepheloid with our turbidity-based method (Section 2.3).
The main shortcoming of the two-basin model is its tendency to overestimate the role of river flushing in removing spilled sediment from Quesnel Lake (Table 3). Only the hypolimnion of the West Basin was highly turbid during the first autumn postspill, so a uniformly mixed model is an especially poor representation of c in the West Basin epilimnion and the Quesnel River during the initial 3 months (Figure 4c). Conditions through autumnal mixis 2014 are somewhat better represented by the model; still, the modeled estimate of cumulative river mass flow through winter 2015, is roughly double the data-based estimate for that period (Table 3), partly because river flow is less than the annual mean (140 m 3 s −1 ), and partly because river sediment c is less than that predicted by the model for the downstream basin (c < c d ). Had the West Basin experienced ice cover during winter 2015, the data-based estimate would have been lower yet. Under ice cover, winter inverse stratification leads to a quiescent surface layer which becomes exceptionally clear as sediment settles out. The surface turbidity sensor on mooring M3 has a nominal depth of 4 m, based on a yearly average. In midwinter, when the water level is at its lowest, this sensor can come to within 2 m of the surface. From mid-February until mid-March 2017, c was near its level of detection (Figure 4a), and the river showed similarly low c (Figure 4c). Based on a particle settling velocity of w s = 5 cm day −1 , clarification would occur to 2 m depth within 40 days in the near-total stillness found under ice. But it depends on the winter: the West Basin may be ice covered for months, in which case ice cover would allow for under-ice clarification; or the West Basin may not be ice covered long enough, or at all, in which case particles suspended in the upper water column would continue to exit the lake via the Quesnel River as they do during the rest of the year.
Apart from seasonal stratification, another aspect of the real system that is not included in the two-basin model is sediment loading beyond the initial spill. Here, we must consider two varieties: autochthonous, originating within the lake due to sediment-water interactions; and allochthonous, originating outside the lake from inflows.
In the case of autochthonous loading, seiche motion has been shown to generate currents that can transport sediment from the bottom into the water column in two ways: through turbulent mixing in the benthic boundary layer and through intrusions (Gloor et al., 1994;Marti & Imberger, 2006;Wain & Rehmann, 2010). Previous studies have included sediment-water interactions in their mass balance models of lake sediment transport (Chapra, 1997;Thomann & Di Toro, 1983). These studies modeled bottom sediments as another CSTR, receiving a flux of sediment from the overlying water column, returning some sediment through the process of resuspension, and losing some sediment to burial. A steady state assumption (dc/dt = 0) allows for the formulation of a net sediment loss rate based on settling, resuspension, and burial velocities. A four-CSTR model could be used to further investigate sediment-water interactions in Quesnel Lake (two basins, each with a bottom sediment layer).
This model would require two resuspension velocity parameters (one for each basin) and two burial velocity parameters, in addition the parameters required for the two-basin model that we have described.
In Section 3.1, we described how snow melt and precipitation during summer 2020 led to increased allochthonous sediment loading. The increase we observed in suspended sediment c in the epilimnion of West Basin we attributed mostly to the Horsefly River. A number of minor inflows also drain directly into the West Basin, and future work on this topic would benefit from a better understanding of allochthonous sediment inputs to Quesnel Lake. During their discussion about potential sources of turbidity in the West Basin, Hamilton et al. (2020) reasoned that Hazeltine Creek was unlikely to be a major contributor of suspended sediment from summer 2015 through 2017, as they had not observed instances of elevated turbidity near its mouth. However, Hazeltine Creek did not flow from its mouth into the West Basin during these years; rather, the creek flowed into a settling pond at its lower reach, then through pipes which discharged into the West Basin at ∼40 m depth via twin, multiport diffusers. From 1 December 2015 until 30 September 2017, MPMC discharged runoff from Hazeltine Creek along with mine effluent through these diffusers; from autumn 2017 onward, the mine piped effluent directly to the diffusers and reestablished the Hazeltine Creek inflow at the shoreline (MPMC, 2018). On 22 May 2018, surface turbidity at ST9 measured over 4 FTU, well above the 0.4 FTU background. The history of mine operation and remediation in the Hazeltine Creek watershed from the spill onward should be assessed in detail when considering allochthonous sediment input to the West Basin in future work.
The usefulness of analytical models lies in their simplicity. The two-basin model provided a broad description of Quesnel Lake's two-phased response to the inflow of suspended sediment following the Mount Polley tailings dam spill, but was not a replacement for data when assessing the actual mass flows of suspended sediment out of the West Basin during the first phase. The analytical approach truly becomes useful over long time periods to describe the net effect of seasonal processes. Assessing the second phase, we found the two-basin model to be most sensitive to settling velocity, so that parameter set the uncertainty in our estimate of 7 ± 5 years to return to background suspended sediment levels. Beyond that time, tailings concentrations in Quesnel Lake will be too low to detect by turbidity alone. More recent observations of turbidity in Quesnel Lake and the Quesnel River support this prediction.

Data Availability Statement
The data used for this analysis are available through Scholars Portal Dataverse (https://doi.org/10.5683/SP2/ RKTIOB).