Seismic Monitoring of a Subarctic River: Seasonal Variations in Hydraulics, Sediment Transport, and Ice Dynamics

High‐latitude rivers are commonly covered by ice for up to one third of the year. Our understanding of the effects of ice on channel morphodynamics and bedload transport is hindered by the difficulties of sensing through the ice and dangers of field work on thin ice or during ice break‐up. To avoid this drawback, we used seismic signals to interpret processes and quantify water and sediment fluxes. Our objective was to determine seasonal differences in hydraulics and bedload sediment transport under ice‐covered versus open‐channel flow conditions using a small seismic network and to provide a first‐order estimation of sediment flux in a Fennoscandian river. Our study reach was on a straight, low‐gradient section of the Sävar River in northern Sweden. Interpretations of seismic signals, from a station 40 m away from the river, and inverted physical models of river stage and bedload flux indicate clear seasonal differences between ice‐covered and open‐channel flow conditions. Diurnal cycles in seismic signals reflecting turbulence and sediment transport are evident directly after ice break‐up. Analysis of seismic signals of ice‐cracking support our visual interpretation of ice break‐up timing and the main ice break‐up mechanism as thermal rather than mechanical. Assuming the bulk of sediment moves during ice break‐up and the snowmelt flood, we calculate a minimum annual sediment flux of 56.2 ± 0.7 t/km2, which drastically reduces the uncertainty from previous estimates (0–50 t/km2) that exclude ice‐covered or ice break‐up periods.


Introduction
Insufficient information on the spatial and temporal variation of the processes under ice-covered conditions limits our understanding of subarctic river dynamics (Ettema, 2002;Lotsari et al., 2019;Turcotte et al., 2011), including estimation of stage-discharge relationships, sediment transport, ice-cover formation, and channel-thalweg alignment (Ettema, 2002;Lotsari et al., 2019;Turcotte et al., 2011). Measurements and thus our understanding of the effects of ice on channel morphodynamics are hindered by the difficulties of "seeing" through the ice and the dangers of making field measurements on thin ice or during ice break-up. High-latitude rivers are commonly covered by ice for up to one third of the year, making investigations of hydraulics and sediment transport challenging. Therefore, case studies from ice-covered rivers have only recently emerged. They have shown that river ice dynamics have the potential to alter cross-sectional flow velocity patterns (Lotsari et al., 2017, shift timing and magnitude of sediment transport (Kämäri et al., 2018), and affect geomorphic change through promoting bank collapse or protecting stream banks during high flows (Tananaev, 2013). Along with seasonal and annual hydroclimatic variations, periodicity of freezing temperatures and the presence of ice cover vary, affecting sediment transport dynamics (Kämäri et al., 2015(Kämäri et al., , 2018Lotsari et al., 2017Lotsari et al., , 2019. These variations control the spatial extent of areas influenced by fluvial versus river-ice processes, which in turn affects the competence of flows, channel bank erosion, and long-term sediment yield. Thus, different seasonal drivers result in different impacts on channel processes and form (Tananaev, 2016). In order to elucidate differences in channel processes, we present a case study of a river in northern Sweden using seismic signals to interpret differences in water flows and sediment flux under ice-covered and open channel flow conditions.

Background
Ice-covered rivers are subject to several geomorphic processes that differ from simple flow-related lateral bank erosion and bed sediment transport. At the beginning of the frozen period, supercooled water can form tiny ice particles (frazil ice) in the water column that attach to the bed to create anchor ice. The presence of anchor ice can significantly increase bedload mobilization thresholds because grain sizes are artificially enlarged through ice-reinforced agglomerate formation, yet allow much coarser sediment to be entrained than through fluvial transport (Kempema & Ettema, 2011;Turcotte et al., 2011). Under ice, flow velocities generally decrease compared to the open-channel flow situation (Ettema & Kempema, 2013;Lotsari et al., 2017), but available data on sediment transport during midwinter stable floating ice cover conditions are limited (Turcotte et al., 2011). However, velocity observations and simulations in a subarctic meandering sandy river suggest that flow characteristics under ice cover are spatially more variable than during open-channel conditions, which should result in greater variability in depositional and erosional locations (Lotsari et al., 2017. River ice break-up combined with the spring snowmelt flood period may have drastic effects on river channel form and sediment transport by erosion beneath ice jams and through ice runs/ice gouging events (Ettema, 2002). Ice break-up can occur in two end-member modes: thermal (a.k.a. overmature) and mechanical (a.k.a. premature or dynamic) (Beltaos, 1997;Deslauriers, 1968). Thermal break-up, through in situ ice decay, may increase sediment load by releasing trapped sediment in ice, and drifting ice pieces can abrade and thus scour stream banks (Turcotte et al., 2011). Mechanical break-up occurs when flow discharge increases before ice decay; more sediment transport is expected in this case due to (1) mobilization of ice frozen to stream banks and the bed, (2) rapidly increasing discharge that generates more turbulence, and (3) thicker ice pieces causing more bank and bed abrasion (Milburn & Prowse, 1996). During ice break-up and winter low-flow periods, plucking of sediment through anchor ice rafting can cause sediment transport of much larger grains (e.g., cobbles weighing 9.5 kg) than is possible based on peak floods (Kempema & Ettema, 2011). Consequently, larger morphologic change (i.e., movement of boulders up to~2 m in diameter) was detected on a gravel/boulder-bed reach of an arctic Finnish river as a result of ice break-up than possible by sediment transport during the snowmelt flood, based on the critical shear stress criterion used for ice-free channels (Lotsari, Wang, et al., 2015). Between these two end-member modes of ice break-up, it is also possible to have ice cracking as a result of increased temperatures; this can also lead to downstream transport of ice blocks with associated sediment, as is common during mechanical break-up. However, researchers commonly consider any fracturing of the ice, even if it has already been subjected to thermal decay, to indicate mechanical breakup (Beltaos, 1997).
Insufficient information on the spatial and temporal variation of the processes under ice-covered conditions limits our understanding of subarctic river dynamics (Ettema, 2002;Lotsari et al., 2019;Turcotte et al., 2011), including estimation of stage-discharge relationships, sediment transport, ice-cover formation, and channel-thalweg alignment (Ettema, 2002;Lotsari et al., 2019;Turcotte et al., 2011). Measurements and thus our understanding of the effects of ice on channel morphodynamics are hindered by the difficulties of sensing through the ice and the dangers of making field measurements on thin ice or during ice break-up. Although ice-covered flow has been studied since the early 1900s (Barrows & Horton, 1907), most existing work related to ice-covered flow and its impacts on hydrodynamics and morphodynamics have been based on laboratory flume experiments (Lau & Krishnappan, 1985;Tsai & Ettema, 1996;Urroz & Ettema, 1994), and technological developments have only recently enabled analysis of ice-covered flows in field conditions. On most rivers, limited methods have rendered the ice season a black box for researchers; thus, we may only see the resulting impacts of river ice but know little about the timing, magnitude, or potential interaction of events. Measuring bedload transport remains elusive, and, compared to suspended sediment (e.g., Kämäri et al., 2018), it is difficult to constrain even in open-channel flow. A typical method for bedload transport measurements use portable traps, such as Helley-Smith or Bunte samplers (Bunte et al., 2004), which require access to the water surface, wadeable conditions in the stream, and for operators to be present in the field, making them difficult to apply in ice-covered rivers. Traditional field methods to measure winter flow velocities, sediment transport, or geomorphic channel change are conducted through holes drilled through the ice and thus only provide point measurements (Demers et al., 2011;Lotsari et al., 2017). Alternatively, before/after comparisons of channel morphology have been used, for example, using high-resolution topography scans (Lotsari, Wang, et al., 2015). None of these prior methods have allowed direct, continuous, high-resolution measurements of flow and sediment transport.
An alternative approach to tackling direct field measurements of sub-ice and ice break-up channel processes can potentially be provided by environmental seismology, an emerging discipline that investigates the seismic signals emitted by processes acting at the surface of the Earth. Previous studies have explored the seismic signals of landslides, rockfalls, debris flows (see reviews by Burtin et al., 2016, andLarose et al., 2015), and also bedload transport and river turbulence (e.g., Burtin et al., 2009Burtin et al., , 2016Roth et al., 2014;Schmandt et al., 2017). These studies have employed portable seismic stations with compact seismometers that are inserted into the ground, connected to data loggers and batteries located above ground. In most studies, fluvial processes, such as turbulent flow and sediment transport, have been interpreted from a predetermined frequency range (e.g., Anthony et al., 2018;Burtin et al., 2008Burtin et al., , 2011Cook et al., 2018;Díaz, 2016), which depends on parameters such as distance to river and seismic wave attenuation properties of the ground. Turbulent flow has been shown to generate a continuous frequency band during perennial flow (~1-20 Hz; Burtin et al., 2008), and sediment transport may appear as discontinuous signals within a higher frequency range (~15-100 Hz; Roth et al., 2016;Schmandt et al., 2017). The exact ranges of these frequencies vary from site to site; therefore, frequency band interpretations are prone to bias, especially when the frequency bands of turbulent flow and bedload overlap. Physical models are now available to predict the seismic spectra emitted by bedload sediment transport and river flow turbulence (Gimbert et al., 2014;Tsai et al., 2012). The models have been tested in both laboratory-and natural-scale experiments and in monitoring efforts by Gimbert et al. (2019) and Schmandt et al. (2017), which found some discrepancies between observations and predicted or anticipated model results, mostly through model underestimation. Proposed explanations for discrepancies include sediment transport mechanisms other than saltation, the application of a model originally designed for bedrock channels to sand and gravel bed systems, and insufficient distance between the source (river) and the seismic sensor. Due to potential for overlapping signals, approximations of solely hydraulics or sediment transport by seismic frequency bands may only work under ideal conditions where the two seismic sources produce minimally overlapping signatures (e.g., Cook et al., 2018). To overcome this ambiguity,  introduced a Monte Carlo-based approach to inversion of the models by Tsai et al. (2012) and Gimbert et al. (2014) for the combined effects of water flow and bedload transport from the record of a seismic station at the bank of an ephemeral stream in Israel. The approach explicitly addresses the spectral overlap by precalculating a reference catalog of combined spectra of potential bedload flux and water level combinations and using this to identify the optimal parameter space to explain the empirical spectra. The authors found reasonable agreement between the hydraulic model predictions and measurements (due to the subordinate seismic power produced by this source) and excellent matches of the model output with independently constrained bedload transport rates, at the uncertainty range of the bedload sampling method itself.
Seismic sensors can also be used to record other environmental sources of seismic waves, such as ice cracking (e.g., Hammer et al., 2015). A small network of seismic stations, installed at safe distances to the active river, can be used to detect and locate discrete cracking events. Therefore, seismic signals may also be used to distinguish between the different modes of ice break-up mechanisms.

Objectives
Here, our main objective was to determine seasonal changes in hydraulics and bedload sediment transport in a subarctic river using a small seismic network, in order to provide a first-order estimation of sediment flux in an ice-covered river in Fennoscandia, the region encompassing the Scandinavian peninsula and Finland. By applying a model inversion approach, we provide a proof of concept for seismic measurements of sediment transport and hydraulics in ice-covered conditions and during ice break-up. Our goal is not to further develop methods as done in previous studies (e.g., Dietze, 2018a; but rather to exploit their predictive capabilities. We compare two different methods of analyzing seismic data to partition signals from turbulent flow and sediment transport. The first method uses a narrow-band proxy, where it is assumed that a continuous signal corresponds to turbulent flow, whereas a higher frequency band corresponds to sediment transport (e.g., Burtin et al., 2008;Roth et al., 2016;Schmandt et al., 2017). The second method uses the seismic model inversion developed by , building on previously developed physically-based models (Gimbert et al., 2014;Tsai et al., 2012). Although, from the onset, we trust the inverse model method over the narrow-band methods, as it is based on physically based processes and not on a somewhat arbitrarily determined frequency range, we show both results to highlight potential differences. The measurement period included mid-winter stable river ice-cover conditions, a transition phase during ice break-up and snowmelt, and open-channel high flow conditions to the end of the falling limb of the snowmelt flood. We make annual sediment flux estimates under the assumption that the bulk of sediment transport occurs during the snowmelt flood, and potentially before and during ice break-up.

Approach and Hypotheses
Within our larger goal of calculating cumulative sediment flux, we aim to address several additional research questions related to hydraulics and sediment transport in ice-covered rivers: (1) How do seismic signatures of hydraulics and bedload sediment transport differ between closed and open-channel flow conditions? (2) How distinct are the average patterns of fluvial processes at high temporal resolution, including possible diurnal cycles, under ice-covered and open-channel flow conditions? (3) Can we detect the seismic signature of the ice break-up and to which extent does it show properties of mechanical versus thermal processes?
We want to test three main hypotheses related to the objectives stated above. First, we hypothesize that there are seasonal and diurnal differences in measured hydraulics and seismic signals of turbulent flow and sediment transport in terms of signal intensity (Gimbert et al., 2014) between ice-covered and open-channel flow conditions (H1). In order to verify differences in flow and potential sediment transport conditions between ice-covered and open-channel flow, measurements are necessary. Acoustic Doppler Current Profilers (ADCPs) enable quantification of the spatial and temporal variation of the hydrodynamics, which contribute to sediment transport, both during open-channel and ice-covered conditions (Demers et al., 2011;Lotsari et al., 2017;Lotsari et al., 2019). However, these measurements cannot show continuous and long-term patterns in hydrodynamics due to the infrastructural constraints of this measurement approach, requiring manual operation. Instead, they can be used as the basis for understanding and interpreting the seismic data sets. During ice-covered conditions, short periods of pressurized flow can occur when gaps in the ice close during freezing; these closed-channel conditions can cause an increase in velocities and localized erosion for a given discharge (Turcotte et al., 2011). Even if flows are not pressurized, flows are confined under the ice with a restricted cross-sectional area, potentially leading to higher velocities and shear stresses. However, in some conditions, sub-ice velocities may decrease due to higher roughness at the water-ice interface (Ettema & Kempema, 2013;Lotsari et al., 2017). Due to any fracturing of the ice (Beltaos, 1997), an ice-covered stream is unlikely completely pressurized directly before ice break-up; therefore, given the same discharge before and after ice break-up, we expect a slight increase in sediment transport and turbulent flow directly after ice break-up. Regarding patterns on shorter spatial scales, diurnal cycles should be evident in turbulent flow and sediment transport seismic signals based on upstream hydrology and snowmelt. These cycles are caused by the diurnal air temperature fluctuations during the snowmelt period (i.e., freezing conditions at night and melting during the day). These diurnal cycles may not be evident until directly after ice break-up if flow is separated above and below the ice cover prior to ice break-up, thus decoupling the seismic signal of turbulent flow from the full flow magnitude above and below the ice.
Second, we hypothesize that bedload transport rates should be temporally synchronous with increases in turbulence, as interpreted from seismic signals, and thus initiate after ice break-up (H2). Because ice break-up could cause a sudden change in velocity and shear stress potential if flow has been separated above and below the ice cover, turbulent flow and sediment transport should increase simultaneously and thus show synchronous seismic signals.
Third, we hypothesize that differences in seasonal timing and intensity of ice-cracking seismic signals can be used to determine the type of ice break-up mechanism, that is, thermal versus mechanical break-up (H3). Ice-cracking seismic signals should appear as nearly instantaneous peaks in frequency and intensity. If mechanical break-up is dominant, then we expect an increase in ice-cracking seismic signals when ice break-up occurs. In contrast, if thermal break-up is dominant, then ice-cracking signals may occur throughout the ice-covered season but should not increase during ice break-up. Although the mechanism behind the formation of ice cracks can be thermal (i.e., melting causes redistribution of stresses), because it is not possible to differentiate between ice crack mechanisms, for our purposes, we assume any ice cracking to indicate mechanical break-up (Beltaos, 1997).

Location
Our study reach was on a straight section of the mostly unregulated Sävar River (a hydropower dam is located~40 km downstream of our site), located 40 km northwest of Umeå in northern Sweden (64.1613°N, 20.3576°E; Figure 1). The location was selected because the reach and upstream section are straight; the cross-section is relatively symmetrical without mid-channel bars, boulders, or islands; and there is minimal instream vegetation that can alter hydraulics during low-flow conditions in mid-winter or summer. Additionally, to carry out sub-ice hydraulic measurements, a thick and stable ice cover was necessary for safety reasons. It was also an ideal site for seismic measurements due to low anthropogenic seismic influence. A one-lane dirt road crossed the river~200 m downstream of the study reach, and a two-lane paved road was located~200 m laterally from the river, on the other side from the seismometers; both roads had minimal and sporadic traffic.
The coastal region of northern Sweden has an average annual temperature of~3°C, with average winter temperatures of approximately −5°C, ranging from −30°C to +10°C, with 500-800 mm of annual precipitation, most of it as snow (SMHI, 2018). During the time of the study deployment until the ice break-up (8 March to 26 April 2018), the average temperature was −3°C (range: −20 to +10).

General Geomorphic and Vegetation Characteristics
The study reach along the Sävar River has a catchment area of approximately 700 km 2 , draining mostly managed forest and some agricultural land. The bankfull width is~20 m, and the maximum bankfull depth is 2.0 m. The channel has a bed slope of 0.0045 m/m and very high relative submergence due to fine bed sediment, causing fairly low-turbulence flow regimes. The grain size distribution on the channel bed is bimodal, consisting of a large proportion of sand/fine gravel and of coarse gravel and small cobbles that armor fine sediment. Based on dry sieving, we calculated a median grain size of 13.5 mm. Winter ice cover contained distinct ice layers: The lower third of the ice was clear and absent of layers of dissolved organic carbon (DOC) or fine sediment, whereas the top half to two thirds of the ice column had an overall light brown color with several distinct bands of darker brown. Based on physical principles, we interpret this to indicate several freezing episodes. The clear portion represents initial freezing, when DOC is expelled from the freezing water (Gold, 1971).
The study reach is situated in a partly confined valley with steeper hillslopes formed within glaciogenic and glacio-fluvial sediment. There is a narrow riparian zone with vegetation (e.g., Carex spp. and Alnus incana) present along the bankfull edge. The upland forest consists of pines (Pinus sylvestris) and spruces (Picea abies) that are able to hold snow on branches that periodically sloughs, especially during wet snow or windy conditions.

General Observations
We set up two time-lapse cameras (one looking upstream and one looking downstream), taking pictures every 15 min. Within sight of the upstream-looking camera, we created a staff gauge on a tree with markings every 5 cm to be able to estimate changes in stage. We also used a calibrated gauge from the Swedish Meteorological and Hydrological Institute (SMHI) (Ytterträsk Nedre, #2236) with 15 min discharge measurements, located~8 km upstream, without any significant tributaries between the station and our study reach. Although it would take water~9 hr to reach our study reach from the gauge, any time correction would carry too many assumptions, regarding diurnal freezing-melting cycles. Therefore, we report the gauge discharge measurements, and we do not think that the use of the daily gauge discharge measurements for calibrating water levels should be affected significantly by a 9 hr lag. We scaled the discharge data (Q w ) from that gauge to river stage at the instrumented site by the Manning equation, where h w is the water level, n is Manning's roughness coefficient, w w is the channel width (20 m), and S 0 is the channel bed slope (0.0045 m/m). Manning's n was estimated through manually minimizing the misfit between model-predicted water level and independently measured water levels. These independent data were either as defined by below the ice cover or based on water levels estimated from pictures of staff gauge during low and moderate water levels. Little research is available on estimating roughness with ice-covered conditions. Therefore, it can be difficult to estimate an accurate n value. As a result, we chose to calibrate n based on independently verified data.
To account for weather events, including wind and precipitation, which can create seismic signals either masking or mimicking fluvial processes, we examined hourly wind and precipitation data from two nearby weather stations: Petisträsk A (64.5669°N, 19.6983°E, 55 km northwest) and Umeå Airport (63.7947°N, 20.2918°E, 40 km south) (SMHI, 2018). Temperature data (every 30 min) was obtained from a logger that we installed at the site.

Topographic and Ice Measurements
During a field campaign with ice-covered winter conditions (6-7 March 2018), we defined a transect across the channel where we measured cross-sectional geometry and ice thickness and conducted measurements of velocity profiles (Figures 1b-1d). Our transect had eleven 20 cm-wide auger holes, ranging from bank to bank and spaced about every 1.5 m. At each drill hole, ice thickness and water depth were measured using a measuring stick and a stadia rod. During open-channel conditions (19 June 2018), we surveyed the transect using a total station (Trimble S8), by wading in the channel, to obtain the channel bed slope, bankfull channel width, and bank geometries.
We used a GoPro camera to take subsurface pictures of the ice layer, which were used to classify the ice undercover roughness (Demers et al., 2013). According to Kämäri et al. (2015), smooth, thermally formed ice cover with a 5-40 cm thickness can decrease average shear stress by~30%, when compared to open water conditions. For rough ice conditions, the differences in shear forces are −1% to + 7% (Kämäri et al., 2015).

Flow and Sediment Transport Conditions
During ice-covered conditions (7 March 2018), measurements of the velocity profile were made in each drill hole using an ADCP, with a SonTek M9 sensor . The ADCP recorded flow velocities throughout the water column in addition to the water depth and channel bed topography. The sensor applies 3.0/1.0 MHz beams, depending on the water depth, yielding an accuracy of ±0.25% (SonTek, 2018). The sensor was mounted on a metal rod and held in place using a tripod ( Figure 1f). A vertical beam was used as the depth reference. Five-minute stationary measurements were performed at each location. Thus, we obtained flow velocity values resolved across the river, vertically from both icecovered and open-channel conditions. The maximum, minimum, and average velocities were obtained from different flow layers in addition to the entire depth-averaged data. We were especially interested in the velocities at the ice-water interface, the near-bed layer, and the depth of maximum velocities. Note that the measurement closest to the ice-water interface is actually 6 cm below. The near-bed velocities refer to the lower boundary of the deepest measurement cell,~2-17 cm above the channel bed, depending on the measurement locations and vertical measurements sampled during 5 min-long measurements. We used a GoPro underwater camera to observe whether there was active bedload sediment transport during ice-covered conditions.
Flow velocities were also measured at the drill hole locations in open-water conditions (18 June 2018) with a moving boat set-up, being held still for 5 min at the hole positions of the winter campaign . In addition, flow velocities and discharge were measured across the section with the moving boat method. The cross-sectional measurements were performed four times to calculate the discharge.

Seismic Instrumentation and Data Processing
We set up a small seismic network ( Figure 1) able to record signals caused by both continuous water flow and sediment transport and by discrete events that can be related to ice cracking activity . Two PE6/B three-component 4.5 Hz geophones were installed in line with the ADCP transect at distances of 10 and 20 m from the right bank (UME01 and UME02, respectively). The signals were logged by Nanometrics Centaur devices recording at 1,000 Hz with a gain of 40. About 40 m upstream of this transect, a Nanometrics Trillium Compact 120 s and a Digos DataCube 3 ext data logger recording at 400 Hz with a gain of 4 were installed~30 m from the right bank (UME04). All sensors were placed 50 cm below the soil-snow interface, and the stations were powered by 12 V AGM batteries. The monitoring experiment lasted from 7 March to 14 June 2018; however, there are some data gaps due to battery or logger malfunctions. For data inversion, we considered data from the most complete data period, 26 March to 7 June.
All recorded data were processed using the R package "eseis Version 0.5.0" (Dietze, 2018a(Dietze, , 2018b. Spectrograms were calculated from the deconvolved and 1-90 Hz bandpass filtered signals using the approach of Welch (1967) with 1 hr-long, nonoverlapping windows, each composed of 360 s-long subwindows with 50% overlap. Following Gimbert et al. (2014), and after inspecting the temporally most stable and outstanding seismic frequency bands (Figure 2a), we used the frequency band of 10-11 Hz as a tentative proxy for river turbulence (i.e., a function of water level assuming all other parameters remaining constant) and the 35-36 Hz band as proxy for bedload transport rate. For the inverse model approach , we calculated sets of 5,000 reference model spectra of the combined effect of water flow (uniformly distributed random water levels between 0.3 and 4 m) and bedload transport (uniformly distributed random bedload fluxes between 0 and 4 kg/s/m), as well as 2,000 reference model spectra of water flow without bedload transport, and combined them to a reference model catalog. These parameter ranges are well beyond the observed water levels ( Figure 3d) and estimated bedload fluxes (Figure 3e). Each of the 1 hr-long empirical spectra were compared to the reference model catalog to identify those spectra with the minimum root mean square differences, along with their deviations (square roots of the squared differences between model spectra and empirical spectra, cf.  from the measured spectra and the corresponding model parameters for water level and bedload flux. We restricted the analysis to the frequency range 5-50 Hz, that is, the part of the spectra that was most significantly affected by river source signals, as is visible in the spectrograms as frequency range with the most energetic continuous signals, and was previously reported by Burtin et al. (2008), Roth et al. (2016), and Schmandt et al. (2017) for other sites (Figure 2b). Since the model output results in strongly fluctuating values, we smoothened them with a running median filter (48 hr filter size) for visualization (Figures 3h, 3i).
Of the 16 seismic model parameters ( Some of the parameters could be independently constrained by field measurements and were considered constant throughout the observation period. These include the average channel width w w (20 m), channel gradient a w (0.0045 m/m), sensor distance to river centerline (r 0 , 20-40 m), and the average grain diameter d s (0.0135 m). The grain size log standard deviation s s (log 0.2 m) was manually adjusted to minimize the root mean square differences of Tsai et al. 's (2012) parametric raised cosine distribution function and the measured grain size distribution of the river bed. Fluid density r w (1,000 kg/m 3 ) and sediment density r s (2,650 kg/m 3 ) were set to typical values for clear water and for quartz, respectively.
Finally, both models require a set of parameters describing seismic ground properties. These include a reference frequency (f 0 , by convention 1 Hz), the ground quality factor (q 0 ), the Rayleigh wave phase velocity (v 0 ), the variation exponent of v 0 with frequency change (p 0 ), and the first of the two Greens function displacement amplitude coefficients (n 0a , cf. Gimbert et al., 2014). Active seismic experiments can in principal be used to constrain these parameters. However, this was not possible due to the snow conditions and an insufficient number of seismic stations available during field work. As an alternative, we modified the parameter combinations manually to minimize the overall deviation of modeled spectra from observed spectra at four different time slices, each as recorded jointly by the three seismic stations. All of the models assume constant seismic ground property values, and the averages used take into account the potential distribution of seismic ground property values, where the natural variability in these values exceed the inter-station variability. This approach limited the parameter space significantly and minimized the potential of equifinality for different parameter combinations. Deviations were quantified using root mean square differences between empirical and model spectra. The four empirical spectra were recorded Journal of Geophysical Research: Earth Surface on (i) 10 March, when there was still a stable ice cover and the water level was determined during field work; (ii) 8 May, during peak discharge as measured by the tree staff gauge (see description in section 3.1); (iii) 28 May, when the measured water level was at 1.5 m; and (iv) 7 June, when the water level was 1.0 m. Fit quality was based not only on minimum deviation but also on how well the individual peaks of the empirical spectra (8-12 and 25-40 Hz) were matched, both in the frequency range and in the seismic power for a given water level. In addition, we required the combined model spectra to resemble the overall shape of the real-world data, in particular the decrease in power toward frequencies <8 Hz and >40 Hz.
The procedure described above allowed us to find a set of optimal model parameters, in terms of minimizing the deviation of model inversion results from independently estimated values. These site parameters are one possible combination of ground properties. It is important to keep in mind that we are not using this model to understand and identify parameter values but using parameter values that best match empirical data. Assuming they adequately represent the average properties of the study system, we used them to invert the seismic data for the target variables.
However, the approach does not account for any uncertainties associated with the optimization of parameters. A further source of uncertainty may arise from site-specific conditions that influence the seismic signals (e.g., Anthony et al., 2018). To account for these effects, we performed the model inversion for all three seismic stations during the time interval when all of the stations had recorded data available. We determined the coefficients of variation (CV; the ratio of the standard deviation to the mean) for both water level and bedload flux during the test period and used the resulting metrics to describe the model uncertainty of station UME04, that is, the station with the most complete record. Note that we expect that the overall model uncertainty is greater than that provided by the multi-station approach, because we were not able to account for the full range of parameter uncertainty.
In principle, the seismic network allows for the detection of discrete seismic activity, such as earthquakes, ice cracks, and road traffic. All of these event types have distinct properties that can be used to discriminate them and isolate the target signals, which in our case were cracks of the river ice cover. Large, teleseismic earthquakes show a sudden onset, distinct arrivals of different wave phases, and a dominant frequency content below 1 Hz. Road traffic, due to its gradual approach and departure, causes a cigar-shaped seismic signal. Both of the signals have durations of several seconds, up to a minute. In contrast, ice cracks, similar to rock cracks (Zeckra, 2015), are sudden, short, impulsive signals with durations shorter than 1 s. To pick discrete events from the continuous stream of data, we applied a STA/LTA algorithm (Allen, 1982), which calculates the ratio of a short-term (0.5 s) to long-term (60 s) running signal average. This method was applied to the envelopes of the 10-40 Hz bandpass filtered signals of station UME04, with an on-ratio of 6 and an off-ratio of 1.5. The waveforms of a randomly drawn subset of 200 detected events were manually screened for plausibility. The sources of the crack signals were located using the signal migration approach (Burtin et al., 2013;Dietze, 2018a) with the same velocity as used for the seismic models described above, applied to a distance grid with a 5 m pixel size. This approach identifies the time delay with which the full waveform envelopes of events are recorded between station pairs and then calculates the cumulative probability density of each pixel in the distance grid for all possible station pairs. This cumulative probability density output was clipped to only include values larger than 95% probability (cf. Dietze et al., 2017), implicitly allowing for source location uncertainty. The clipped location estimates of all events were normalized and summed, which yielded a spatially resolved estimate of overall ice-crack activity. We use this location map predominantly to ensure that the locations of the signals are not distributed randomly or concentrated along the road as a potential anthropogenic source.  (Gimbert et al., 2014) and Bedload Sediment Transport (Tsai et al., 2012) Parameter (

10.1029/2019JF005333
Journal of Geophysical Research: Earth Surface

Sediment Flux Calculations
We estimated hourly sediment fluxes using the model inversion method described above with Tsai et al.'s (2012) model. The average specific sediment flux (m 3 /sm) was converted to total sediment load for each hour, standardized by the catchment area (t/km 2 ), using a channel width of 20 m, sediment density of 2,650 kg/m 3 , and drainage area of 700 km 2 . Finally, the hourly sediment loads were summed for the period of record to estimate the cumulative sediment flux. Assuming a sediment density of 1,700 kg/m 3 , taking porosity in to account, we transformed the cumulative sediment flux to an estimate of the catchment-wide erosion rate. The estimated error on the cumulative sediment flux was calculated using Gaussian error propagation for each hourly value, according to Equation 2. The error value was calculated as the median CV value based on the seismic stations that provided reliable data.
where EE is the estimated absolute error, CV is the coefficient of variation, and q s is the hourly sediment flux in t/km 2 . In addition, we estimated the order of magnitude of expected bedload flux using the modified Meyer-Peter Müller equation (Equation 3; Wong & Parker, 2006) based on flow depth from the scaled discharge data.
where q s is the specific bedload flux (kg/s/m), d s is the median grain size, h w is the water depth, S 0 is the channel slope (0.0045 m/m), r s is the sediment density (2,650 kg/m 3 ), r w is the water density (1,000 kg/ m 3 ), and θ c is the critical dimensionless shear stress (0.0495). The Meyer-Peter and Müller (1948) equation has not been validated for ice-covered streams, and we do not expect it to give a reliable prediction. It is used here to cross-check the order of magnitude obtained from the seismic inversion.

Hydroclimatic Observations
During the study period in 2018, temperatures ranged from −20°C to +10°C prior to ice break-up, and from 0°C to 30°C after ice break-up; wind speeds typically ranged from 0-3 m/s with a peak at~5 m/s during a wind storm event in early June (Figures 3a-3c). During the winter field survey (7 March 2018), there was a stable, snow-veneered ice cover on the reach. The ice roughness of the ice-water interface was visually interpreted as "smooth" ice with gentle undulations, following the definitions by Demers et al. (2013). This ice surface smoothness indicates low ice-water interface roughness, and thus lower roughness values than the channel bed with gravel and cobbles. There was no frazil or anchor ice observed (based on go-pro imagery and visual observation through drilled holes) during the field survey. During mid-winter stable ice cover conditions, the ice thickness ranged from 0.41 m (Hole 6) to 0.65 m (Hole 9), with an average thickness of 0.56 m (7 March 2018). Approximately 0.15-0.20 m of fresh snow resided on top of the ice.

Ice Break-Up Timing and Types
Static ice conditions lasted until mid-April, based on time-lapse imagery and maintenance field visits ( Figure S1). On 6 April 2018, a steeper, more turbulent reach,~150 m downstream of the study site, showed patches of open water. The study reach still had a stable ice cover and was completely closed with no visual changes in ice conditions. During a 21 April 2018 site visit, the downstream turbulent reach was mostly open, but stream banks still had ice and snow cover in addition to ice bridges. The open section had extended about 20 m upstream toward the study reach. The study transect had a thin and wet ice cover in the middle of the channel, indicating that thermal melting had begun and that the flowing water had likely excavated the ice cover from the bottom of the surface ice. On 22 April 2018, about 10 cm of water was flowing on top of the ice, most likely originating from open water upstream. This indicates flow separation above and below the ice and that the sub-ice flow conditions were not pressurized; thus, the entire discharge was not able to exert shear forces on the river bed. Thermal melting continued both due to the high air temperatures and flow

10.1029/2019JF005333
Journal of Geophysical Research: Earth Surface eroding the ice from the top and bottom of the ice cover. On 25 April, thermal decay progressed as the studied cross section was progressively opened along a patch of open water extending downstream and prograding sideways. Thermal versus mechanical break-up was identified in photographs as in situ melting, where the open area extends laterally, or as distinct cracking and movement of ice blocks, respectively. At noon on 26 April, the ice frozen to the banks had fragmented and started floating downstream. On 27 April, the entire section was free of ice. During the end of April, there was still a consistent snow cover in the surrounding forest, but it was wet and shallower than during the March field campaign. Snowmelt water had percolated through the snow into the ground, and melting was only visible around the tree trunks.

Seasonal Hydrodynamics
The discharge during the ice-covered winter measurements ( When the 5-min measurements were time-averaged (as measured at each drill hole location), the ice-water interface and depth-averaged velocities were lower in ice-covered conditions than during open-channel flow conditions. During winter, the high-velocity core was located in the middle of the water column rather than at the surface as seen in open-channel flow conditions. This maximum time-averaged velocity was~9 m from the left bank (36.6 cm/s, Hole 5; Figure 4a). In winter, the depth-averaged (and time-averaged) velocities of the whole water column ranged from 8.9-32.9 cm/s (not shown in Figure 4a), which was close to the range of time-averaged minimum and maximum velocities of the ice-water interface (7.6-32.1 cm/s; Figure 4) and the near-bed layer (7.1-29.6 cm/s; Figure 4a). The highest winter depth-averaged velocity was at Hole 5 (32.9 cm/s), located in the middle of the channel, whereas the lowest depth-averaged flow velocity occurred along the right bank at Hole 11 (8.9 cm/s). The greatest maximum time-averaged velocity was in the middle of the water columñ 9 m from the left bank (36.6. cm/s, Hole 5; Figure 4a). The greatest ice-water interface velocity was located at Hole 6 (28.5 cm/s; Figure 4a). The lowest near-bed, maximum, and ice-water interface velocities were at Hole 11 at 7.1, 10.9, and 7.6 cm/s, respectively.
In contrast, in the summer (June 2018; Figure 4b), the highest depth-averaged (and time-averaged) velocity (Hole 8, 28.1 cm/s; note that the depth-averaged values are not shown in Figure 4), near-bed layer velocity (24.6 cm/s, Hole 10; Figure 4b) and surface velocities (31.6 cm/s, Hole 8; Figure 4b) occurred closer to the right bank than in the ice-covered flow conditions. The high-velocity core in open-channel flow was closer to the surface than when ice cover was present (Figure 4b). However, the largest of the maximum time-averaged velocities (32 cm/s, Hole 8; Figure 4b) was slightly less than that in the winter. In most of the measurement columns, the surface layer velocities were greater than the depth-averaged velocities. The minimum depth-averaged (and time-averaged) velocity was 7.7 cm/s, located next to the left bank. The minimum surface velocity and the minimum near-bed velocities were 8.7 and 3.7 cm/s, respectively, which were all located next to the left bank, that is, measured at Hole 1 (Figure 4b).
Based on the Go Pro camera videos of the channel bed recorded on 7 March 2018, there was minimal to no sediment transport. No sediment was in motion on the bed close to the banks, which consisted of pebbles and gravel. Only a few separate grains were in noncontinuous motion in the middle of the channel, but the bed was likely mostly static under the ice-covered midwinter conditions. Unfortunately, we were not able to conduct similar observations during the summer campaign due to an equipment malfunction.

Continuous Seismic Measurements
The seismic data are mainly interpreted from their time-spectral perspective. We focus our report on data from station UME04 because it comprises the full study period without any data gaps. The spectrogram (Figure 2a) of the broadband station UME04 (see Figure S2 for plots of the other two stations) shows three continuously active frequency bands. The frequency band below 2 Hz is excited predominantly during periods of high wind speed. Another frequency band between 8 and 12 Hz shows a diurnal pattern of energy state, suppressed before and intensified after ice break-up. Superposed on the 8-12 Hz band is an increased seismic energy between the end of April and beginning of June. A third frequency band exists at 20-50 Hz, which does not emerge until late April, after the 8-12 Hz band had already shown increased energy. From 3 May onward, there is a strong diurnal pattern present in this frequency band (visible as daily humps around 40 Hz; Figure 2a inset), progressively decreasing in intensity until the pattern is no longer discernible by the end of May. From time to time, broadband bursts of seismic energy appear (vertical patterns in Figure 2a), which usually coincide with windy days (Figures 3b, S3b, and S3h). The most prominent of those events is on 4 June (Triangle 3 in Figure 2a), when the wind speed record shows its maximum values. Two additional events, with slightly narrower frequency bands with energy mostly between 8 and 35 Hz, are also visible in the spectrogram: on 20 April and 25-26 April (Triangles 1 and 2 in Figure 2a).
The tentative proxies for discharge (10-11 Hz band; Figure 3f) and bedload transport (35-36 Hz band; Figure 3g) appear to be decoupled from the evolution of imagery-based and modeled river stage (Figure 3d). Instead of showing the linearly increasing trend from 26 April to 8 May and the slowly falling limb of the flood for the rest of the time, the 10-11 Hz band rises and falls more gradually with a convex shape that reaches a maximum around 20 May. The 35-36 Hz band indicates only subtle changes over the entire time, after its state shift during the ice break-up phase.

Seismic Model Inversions
The best fit combinations (i.e., parameters yielding the smallest deviations and spectral shapes of combined water level and bedload models for multiple seismic stations and time slices of environmental conditions) for the unknown parameters were ground quality factor q 0 = 11 and Rayleigh wave phase velocity v 0 = 700 m/s. These results are in the range of expected values (e.g., Boxberger et al., 2017) for the sand-to cobble-till deposits and finer alluvial sediment present in the study area. The variation exponent of Rayleigh wave velocities p 0 was 0.79, and Greens function displacement amplitude coefficients, n 0a = 0.6 and n 0b = 0.8 (Figure 2b). In addition, the exponent characterizing the quality factor increase with frequency, e 0 , had to be set to 0.49; otherwise, the model spectra did not converge to a meaningful fit (for a discussion on equifinality of parameterization, see section 5.1). Figure 2b shows how the above parameter combination principally resembles the frequency modes and tails of the empirical spectra of the four different test water levels, despite having an overall much smoother and wider distribution of the frequencies (see Figure S3i for the corresponding fit errors). The figure also shows that high-quality factors (q 0 ) may lead to shifts in the frequency peak position away from the empirical ones, overestimation of the seismic power, and a lack of decreasing power for frequencies >40 Hz. Likewise, too low Rayleigh phase wave velocities resulted in drastic shifts of the peak frequencies and amplitudes away from the empirical data.
Until the emergence of the 20-50 Hz band (Figures 2a and S3h), the model inversion fit error ( Figure S3i) is high for frequencies around 10-15 Hz and below 5 Hz. After that emergence, misfits mainly appear above 25 Hz. All three seismic stations produced overall comparable spectral data, with systematically lower seismic powers and lower active frequency bands with greater distances from the stream (Figures 5a and S2). At a first glance, inversions of water level from all three stations were overlapping, with average CV of 0.25 (Figures 5b and 5d). However, upon closer investigation, the model inversion results for Station UME01 show that the inverted values basically only yielded discrete water level values of 0.6, 1.3, 2.0, and 2.3 m with almost no solutions in between these clusters. Likewise, the transition from low winter flow toward increased spring flood levels is not represented by the UME01 water level model inversions (black curve in Figure 5b, left plot). Inversions of bedload fluxes were similar between UME02 and UME04, whereas UME01 values were almost exclusively zero, with only short excursions toward values an order of magnitude below the values of the other two stations (Figure 5b, right plot). The expected increase in bedload flux values during the spring flood was not visible in the modeled values from station UME01. Since these stark and outstanding differences of this single station cannot be explained in a straightforward way, we excluded 10.1029/2019JF005333 Journal of Geophysical Research: Earth Surface this data set from the calculation of the CV. The calculated CV values, calculated from UME02 and UME04, were 0.25 for water level and 0.38 for bedload flux. Because the distribution of CV values was highly skewed, we used the median rather than the mean. Since data from station UME01 were rejected, due to unrealistic model results, and station UME02 only recorded until the end of April, we based our model inversions on Station UME04, which provided continuous seismic data. The CV values were used to provide a first-order estimate of the model result uncertainties (Figures 3h and 3i). After smoothing the raw model inversion time series with a running median filter (48 hr filter size), we added the ±25% uncertainty range to the water level and ±38% uncertainty range to the bedload flux time series (Figures 3h and 3i).
The estimated average water level time series obtained from the inverse model is in agreement with the time-lapse image-based and Manning-scaled determinations of water level (Figures 3d and 3j, R 2 = 0.86 and 0.93, p < 0.0001). Manning's n, which was calibrated based on independently obtained data (Q w and

10.1029/2019JF005333
h w ), was estimated at 0.2. Although this is quite high for this channel type (Chow, 1959, suggests 0.03-0.04), the value allowed us to relate the discharge from the upstream gauge to water levels measured at our site; lower values yielded lower than realistic values of water level based on independent data from time-lapse photos of the staff gauge and the measured water level below the ice. The Manning-based model overestimated the average depth of the freely flowing water column during the winter field campaign (0.8 m) by about 29 cm on average, and the peak water level during the flood (2.1 m) was underestimated by 0.4 m compared to the camera-based maximum value (2.54 m). During 25-26 April, the inverted water level shows a spike, coinciding with a burst of seismic energy in the 8-35 Hz band (Triangle 2 in Figure 2a).

Sediment Flux Calculations
Modeled bedload fluxes are between 0.01 and 0.04 kg/s/m during the ice-covered period, start rising from around 22 April and range between 0.5 and 0.7 kg/s/m between May and early June (Figure 3i), until beginning to recede toward the end of the recording period. A spike occurs around 5 June (vertical bar in Figure 5). The range of seismically modeled bedload flux is at the same order of magnitude as predicted by the simple Meyer-Peter Müller estimate (Figure 3e). However, the temporal evolution differs significantly. Based on the seismically modeled sediment flux (Figure 3l), the cumulative sediment export for the period of record (26 March to 7 June 2018) was 59.7 ± 0.7 t/km 2 /a, of which 2.8% occurred under ice-covered conditions (before 22 April), 1.7% occurred during the ice break-up (22-26 April), and 95.5% occurred during or after the snowmelt flood. Assuming a bulk density of 1,700 kg/m 3 , this sediment flux translates to a catchment-wide erosion rate of 35.1 ± 0.4 μm/a.

Ice Crack Analysis
The STA/LTA algorithm yielded a total of 4,532 potential events. Their duration distribution peaks at 0.5 s. Around 0.2 and 2 s, the distribution has fallen off to 10% of its maximum height ( Figure S4). Thus, we rejected all picks that lasted less than 0.2 or more than 2 s. Inspection of the randomly drawn subset (sample size of 200) of the remaining 1,553 events revealed that almost all events showed a sudden onset of the signal, followed by a gradual decrease toward background noise. Most of the energy is contained in the 10-90 Hz frequency band.
The source locations of these crack signals are predominantly clustered along the river (n = 899, Figure 6a) with some isolated smaller clusters around the bridge, 150 m downstream of the transect, and further northwest in the forest. Although the network geometry is not ideal for locating seismic sources, the results allow us to rule out the road as the dominant source of the identified seismic signals. Crack activity by season (Figure 6b) shows a major crack period on 26 March, during which there were large temperature variations (−2 to −18°C). Between the beginning and end of April, the average crack rate was 45 events per day. The crack rate suddenly decreased to about 5-7 events per day after 26 April. Crack activity is highest in the early afternoon (Figure 6c).

Seismic Signals Reveal Seasonal Turbulent Flow and Sediment Dynamics
We used seismic signals to generate a first estimation of sediment flux in an ice-covered subarctic Fennoscandian river, which provides both a geographical area and a geomorphic process lacking sediment flux data. In high-latitude rivers, seasonal changes in turbulent flow and sediment transport can be affected by hydrological changes due to snowmelt-flood peaks as well as changes in ice dynamics. Due to difficulties

10.1029/2019JF005333
Journal of Geophysical Research: Earth Surface measuring these processes around the time of snowmelt flood peaks and ice break-up, past studies have been unable to parse out the effects of hydrology compared to ice dynamics and removal of ice cover. Based on velocity data, and interpretations of seismic signals and their modeled results, we found clear seasonal and diurnal differences between ice-covered and open-channel flow conditions, supporting H1 that there are seasonal differences in hydraulics and seismic signal properties due to ice-covered versus open-channel flow conditions. Winter velocity profiles exhibit a high-velocity core, based on time-average data, in the middle of the water column and toward the left side of the channel (Figure 2). However, the highest instantaneous velocities were located toward the right side of the channel, indicating very large spatial and temporal heterogeneity in velocities under ice-covered conditions. The open-channel flow conditions exhibited the more typical high-velocity core at the water surface, located toward the right side of the channel. This sudden shift in the location of the high-velocity core may also be evident in seismic data as a stepped increase to a slightly broader frequency of turbulent flow signals. However, we recommend more research that focuses explicitly on the differences in ice-covered and open-channel flow under similar flow conditions. We interpret that the seismic spectrogram of station UME04 (Figure 2a, like spectrograms from other stations) is dominated by fluvial sources. While the 8-12 Hz frequency band persisted throughout the entire study period, that is, also during ice-covered conditions, the 20-50 Hz band only exhibited seismic power with the onset of the spring flood and thus had an increased likelihood of sediment movement at the river bed. With the receding spring flood in June, the 20-50 Hz frequency band started to diminish again, while the 8-12 Hz band persisted. This clearly indicates two different seismic sources (see also similar findings of time evolution and frequency bands by Schmandt et al., 2017), with their power driven by changes in discharge during this time of the year (i.e., visible as water level in our independent data). Exploration of the model space (Figure 2b) illustrates the sensitivity of seismic ground properties such as quality factor and Rayleigh wave phase velocity. Thus, it is crucial to estimate these parameters with sufficient rigor, based on several time slices with different system states (i.e., water levels and assumed activity/absence of bedload transport) and multiple seismic data sets. A value of 0.79 for the variation coefficient of the Rayleigh phase wave velocity (p 0 ) imposes a strong decrease of wave velocities for frequencies above about 40-50 Hz. However, only with this model setup were we able to find reasonable fits with the empirical data. Thus, in this study area setting, the bedload model is either confronted with different seismic source processes responsible for bedload transport or with a different wave propagation mechanism. In any case, with the instrumentation we used, this cannot be investigated with enough detail to finally resolve the question. There is an indication from other studies that there are some limitations to the seismic model (Gimbert et al., 2019). However, in this study, our aim is to interpret Earth surface processes by applying the most state-of-the-art models, but we encourage other researchers to continue polishing the mechanics of the models.
In the absence of active seismic survey data, our approach was to minimize the deviation between empirical and modeled spectra for characteristic time slices, based on the combined results of multiple seismic data sets. This yielded an overall range for the target variables water level and bedload flux of 25% and 38%, respectively. These values are higher than the results shown by  and suggest that more than two to three sensors should be used in future studies to define the uncertainty range of model outputs (as shown in Figure 5). But even then, since our data did not cover the entire observation period with all stations, our uncertainty estimates are limited. Further uncertainties are introduced by using fixed values for model parameters such as q 0 , v 0 , or p 0 . The inverse model approach would in principle allow for extended Monte Carlo lookup tables, which also vary parameters other than water level and bedload flux. Thus, while in principal there may be other parameter combinations that yield similar agreements between combined model results and independent data (equifinality), it is beyond the scope and data in this study to explore such phenomena in detail. Although the combination of parameters is one out of potentially more solutions, we optimized for a solution that describes the reference water levels and the overall shape of the empirical spectra with the least deviations. An active seismic survey is necessary to adequately constrain these elements of the model. As a consequence, the uncertainty estimates reported here (25% and 38% for water level and bedload flux, respectively) need to be considered as minimum estimates, only accounting for interstation variability.
During the uncertainty analysis, we decided to reject the records from Station UME01. This decision was based on the obviously spurious results of the model inversion (Figure 5b), which did not agree with any physically grounded expectations about the evolution of water level and bedload flux (see section 4.5 for further description). Furthermore, the gradual increase in water level after the ice break-up was not depicted by the inversion data, and the inverted bedload flux time series did not show any systematic commonalities with the other two stations nor, more importantly, with the rising water levels that would trigger an increased flux of bedload particles. Other studies have discussed the susceptibility of high-frequency signals to local effects (e.g., Anthony et al., 2018).  showed that at small distances to the river, spatially heterogeneous sediment flux can bias the model output. In addition, the spectral data of Station UME01 (Figure 5a) do not show such a clear separation of the two seismic sources with respect to the frequency ranges covered, as seen in UME02 and UME04, due to the different attenuation relationships of the two sources with distance. Thus, we are not surprised by the finding that the station that is closest to the seismic source of interest was most affected by these above-mentioned effects. The other stations, located two and four times further away, were affected by these issues to a smaller degree, and thus, their records provided more homogeneous model results. Thus, installations a medium distance from the channel may be able to more accurately model water levels and sediment flux. Despite the structural limitations of the seismic models (cf. discussions by Anthony et al., 2018;Gimbert et al., 2019;Schmandt et al., 2017) and the observed deviations, the Monte Carlo-based inversion approach linking both models was able to return meaningful results that are in agreement with independent data on water level and calculations of sediment fluxes (Wong & Parker, 2006; Figure 5).
On average, the model overestimated water level during the ice-covered period. A reason for this effect may be the confined flow conditions, leading to a different distribution of vertical flow velocities and thus higher than expected seismic energy transmission into the ground. Likewise, the zone of highest flow velocity was closer to the seismic sensors and closer to the bed in winter compared to the summer survey (Figure 4), leading to more seismic power recorded for the same average water level, simply due to changes in the source-receiver distance (see also discussion by . However, the uncertainty range of the model output is too high to systematically explore reasons of this effect with sufficient rigor, and we leave this point a topic of future experiments. The model underestimates water level for high discharges. We interpret this as the result of a step change in fluvial dynamics when the water level rises above bankfull depth (about 2 m). Under these conditions, the river width, grain-size distribution, flow velocity, and average depth change dramatically. The model does not account for these changes. Similar trends are visible when comparing the Manning-scaled water levels with the image-based determined water levels. The latter also yielded lower values when the river rises above 2 m (Figure 3d). Comparison of the modeled bedload flux (Figure 3i) with the Meyer-Peter and Müller (1948) approach (Equation 2) yields broad agreement, although the latter approach only serves as an order of magnitude estimate rather than a precise proxy of temporal changes in bedload flux under natural conditions. Thus, here we do not observe an order of magnitude difference between modeled and observed seismic power due to bed particle impacts as has been reported by, for example, Schmandt et al. (2017). This difference may be explained by the differences in channel configuration at our site at the Sävar River, and the Trinity River reach instrumented by Schmandt et al. (2017). While the Sävar River has an equilibrated straight channel with a sand and gravel bed, the Trinity River features a much more active gravel bed channel, into which material was injected from a discrete source proximal to the seismic sensors, which was mobilized by a dam-controlled flood for just a few hours. This operation imposed transient conditions to the river system, resulting in local rearrangements of the system, for example, by spatially heterogeneous bedload flux and bar migration.
The simple approach of using the 10-11 Hz frequency band (Figure 3f) or 35-36 Hz band (Figure 3g) to describe the water level and sediment flux, respectively, does not give meaningful results in comparison with any of the other utilized measurement or model data. Whereas the hydrograph clearly depicts the spring flood and its receding limb, the 10-11 Hz frequency band continued to remain high with only a slight receding trend. Likewise, the 35-36 Hz band (Figure 3g) did not provide any interpretable process-driven insight, as it did not show any relationships with water level. Thus, the temporal evolution of both proxy time series does not show agreement with the documented fluvial conditions. The ice break-up coincided with an increase in discharge magnitude, as is common in many high-latitude rivers, but the intensity of seismically sensed turbulent flow and sediment transport also reflects changes in ice conditions. There is a sudden increase in seismic power at the time that ice cover melts and completely disappears (Figures 3h and 3i; see gray line depicting raw model output), which may either be due to a sudden, short-term increase in turbulent flow and sediment flux or potentially an artifact of increased levels of ice cracking. An increase in turbulent flow and sediment transport after ice break-up would run contrary to what would be predicted based on physical models and processes (Ettema, 2002;Lotsari et al., 2017). Changes in hydraulic conditions after ice break-up are, according to such models, caused by either transitions from closed-channel flow (e.g., pipe flow) conditions, with higher pressures and thus shear stress, to open-channel flow, or due to ice cover causing confining flows and thus increasing velocities and shear stress. However, Turcotte et al. (2011) have found that even minor cracks in the ice cover cause nonpressurized flow conditions under the ice cover. Through the time-lapse photos, we observed ice cracks either within our transect or upstream and downstream, suggesting that the Sävar River had not experienced pressurized flow conditions under ice at the study site. Another potential mechanism for increased turbulent flow and sediment transport after ice break-up is separation of flows above and below ice prior to breakup. During the rising limb of the hydrograph, flow was observed both below and above the continuous ice cover, which could have originated from upstream ice break-up or local cracks. If the flow is separated above and below the ice, then the shear stress and velocities, and thus seismic signals, will not reflect a commensurate increase in discharge.
Diurnal cycles in seismic signals reflecting turbulent flow and sediment transport are evident starting directly after ice break-up and increase in intensity with increased flow magnitudes (Figure 2a). The diurnal pattern persists throughout the rising limb and most of the receding limb of the snowmelt flood. This pattern in the seismic data, in particular within frequencies related to sediment transport, mirrors patterns seen in 15-min flow data from the upstream gauge. Between the time of ice break-up and the peak flow, there is an inverse relationship between flows at the upstream gauge and calculated sediment flux ( Figure S5). Such high-resolution temporal data of bedload transport is nearly impossible to obtain using traditional field methods.

Potential Asynchronicity in Sediment Dynamics in Ice-Covered Rivers
Traditional methods for determining the onset of sediment transport and quantifying its magnitude have not allowed measurements during or directly prior to ice break-up. Using seismic measurements, we have been able to test Hypothesis H2 that sediment transport seismic signals should be temporally synchronous with seismic signals of turbulent flow and ice break-up. Based on basic interpretation of narrow frequency bands (i.e., 10-11 Hz for turbulent flow and 35-36 Hz for sediment transport) from the seismogram, it appears as if sediment transport signals were asynchronous with that of the seismic signals of turbulent flow and ice break-up. In contrast, the bedload transport calculations and river stage from the inversion models show parallel evolution of these two processes when examining raw model output (Figures 3h and 3i), including a significant increase in both bedload flux and turbulent flow occurred just prior to ice break-up. Because of overlapping effects of turbulent flow and bedload signals, the use of simple frequency bands (e.g., 10-11 Hz or 35-36 Hz) as proxies to interpret the onset or strength of hydraulic parameters may be misleading. Therefore, fine-scale patterns in turbulent flow and bedload seismic signals based purely on frequency bands may be strongly influenced by bed roughness, which may change over the course of a flood (Roth et al., 2017). Thus, we recommend using physical model inversions  to interpret fluvial processes of turbulent flow and sediment transport.

Seismic Interpretation of Ice Break-Up Mechanism
Analysis of seismic signals of ice cracking support our visual interpretation for the timing of ice break-up and of the main ice break-up mechanism as thermal rather than mechanical, supporting H3 that we can predict break-up mode using the temporal distribution of ice cracking during the winter and ice break-up period. In this analysis, we define all ice cracking as indicative of mechanical ice break-up. The cause of ice cracking, either by mechanical forces or temperature changes, is not taken into consideration. However, the diurnal distribution of cracking (Figure 6c), with most cracks occurring in the early afternoon, suggests that most cracks are thermally induced. Ice-crack signals were seen in the seismic data throughout the ice-covered measurement period (March to April), and they were generally located within or very close to the channel, confirming that ice cracking is the most likely source process (Figure 6a). The few isolated clusters away from the river might be artifacts of the location approach or particularly strong, short impulsive events by sources other than ice cracking. Although we cannot resolve the role of these secondary small hot spots, the main activity cluster matches locations close to or on the river. The temporal distribution of ice-crack signals was fairly even throughout the ice-covered period, except for a large sudden increase at the end of March, on a day with large temperature variations (Figure 6b). If mechanical break-up was the main mechanism, we would expect to either only detect ice-crack signals during or right before ice break-up or detect a sudden increase in the number or intensity of ice-cracking signals at this time. No large increase in the number of ice cracks per day occurred at the time of ice break-up, except for the last day, during which ice blocks were present and mechanical break-up was evident in time-lapse photos (Figure 6b). Thus, the ice-cracking signals are consistent with the image-based interpretation of ice break-up on the Sävar River in 2018, which started as thermal break-up, showing no increase in ice-crack signals. A minor mechanical element is reflected by the increasing number of ice-cracking signals during the final stages of break-up. Furthermore, ice-cracking signals decreased drastically in number after all ice had disappeared from the reach, grading into a background rate due to the noisy environment. This provides an independent method of determining the timing of ice break-up if photographic evidence is lacking.
The mode of break-up has implications for geomorphic and ecological processes (e.g., Lind et al., 2014). For example, mechanical break-up potentially leads to a larger increase in sediment transport intensity that should be synchronous with ice break-up (Milburn & Prowse, 1996). Therefore, having a method to easily determine the degree to which the ice break-up is steered by these end-member processes is vital for understanding river ice dynamics. Currently, other methods for determining ice break-up methods are visual observation and using a shallow water ice profiler. The latter gives point measurements of ice thickness but involves extensive instream installation efforts and costs. Visual observations, which are labor intensive if made in person or by sorting through photographs, can be unreliable for several reasons: (1) Observation gaps between photographs or during periods of low visibility (e.g., fog, night, and snow storm) can omit important processes, and (2) snow cover can cause misleading interpretations, as small ice cracks that can lead to mechanical break-up may not be visible. In contrast, seismic sensors can be installed along the banks of a larger reach to track the break-up processes as they propagate downstream.

Sediment Flux of Subarctic Ice-Covered River
Using continuous seismic measurements, we could calculate a first-order estimate of sediment flux during the most dynamic time period of a subarctic river. Because the Sävar River has a snowmelt-dominated flow regime, with some lower-magnitude flashy rainfall-induced discharge peaks during autumn, the snowmelt flood should transport the majority of sediment. By applying the 48 hr running median filter to the data (bold line in Figure 3i), we can account for and remove short periods of increased seismic activity that are coincident with meteorological events. Thus, only the multiple day event on 3-6 June affects the results in a significant way. When removing the data for this period and interpolating the bounding data, the result changes from 59.7 ± 0.7 t/km 2 /a to 56.2 ± 0.7 t/km 2 /a. Thus, assuming that the bulk of sediment moves during ice break-up and the snowmelt flood and is transported as bedload, our result of 56.2 ± 0.7 t/km 2 /a can be considered as an annual estimate. In any case, it gives a minimum estimate for sediment flux in this subarctic Fennoscandian catchment controlled by a snowmelt flow regime and ice cover and break-up. As such, it is one of the first reported sediment flux values of the most dynamic time of the year, taking into account the full rising limb of the snowmelt flood during sub-ice and ice break-up conditions. Our values generally agree with values calculated from the Meyer-Peter and Müller (1948) sediment transport equation (Equation 2). Therefore, seismic measurements used to model sediment flux provide an alternative method for determining sediment transport and hydraulics in ice-covered conditions and during ice break-up.
Based on our modeled sediment flux, the catchment-scale erosion rate corresponding to the sediment export of 56.2 ± 0.7 t/km 2 /a at our site is 35.1 ± 0.4 μm/a, which is 1 order of magnitude lower than in cold mountainous regions (e.g., the European Alps: 0.3-1 mm/a; Delunel et al., 2010). Very few reported values of sediment flux in similar catchments exist, and estimates from global reviews of sediment flux consistently place Fennoscandia in the lowest bin of values: 0-50 t/km 2 /a (Walling & Webb, 1983) or 0-5 t/km 2 /a (Lvovich et al., 1991). These low values reflect a landscape with crystalline bedrock, relatively low relief, and a cold 10.1029/2019JF005333 Journal of Geophysical Research: Earth Surface climate contributing to low rates of weathering and erosion. Since field data have been lacking from ice-covered periods, previously reported values may have underestimated annual sediment flux. In the Sävar River, nearly 5% of the measured sediment flux occurred before open-channel flow conditions-a period that was previously unmeasured.

Conclusions
Our results show strong potential for using seismic instruments for understanding fluvial geomorphic processes during the previously understudied winter season at northern latitudes due to the difficulties of field measurements under ice or during ice break-up. We have shown a proof-of-concept for estimating sediment flux. Fine tuning with field calibration of both seismic ground parameters and sub-ice sediment transport may be required to answer questions of, for example, the role of ice dynamics in shaping channel form and changes in sediment transport during thermal or mechanical ice break-up. Future studies will need to calibrate sub-ice bedload sediment transport with seismic data. Although methodologically difficult, some options exist to obtain rough estimates of bedload transport, such as Helley-Smith samplers, which can be used in drill holes during stable ice conditions . Either an in-stream acoustic system or fix-installed slot samplers can be used, but they have both a significant impact on the system and require high installation and maintenance efforts. Furthermore, bedload data should be collected over a range of conditions to capture variations in ice seasons with different types of ice break-up, ranging from completely mechanical to completely thermal. Seismic data of the ice season from more geomorphic settings, including a range in grain sizes, bank types, and channel slopes, would allow us to constrain in which environments ice processes play a primary role in sediment transport and thus determine channel form. The timing of bed sediment transport and bank erosion in relation to ice break-up relative to the snowmelt flood has in practice been impossible to constrain due to difficult and life-threatening conditions for fieldwork. But with the remote and continuous nature of seismic sensors, we should be able to parse out the separate roles of ice jams, ice break-up, and flow-induced sediment transport and erosion in the total sediment yield. Therefore, we strongly recommend the use of seismic methods in ice-covered streams; however, some difficulties can arise when there is abundant anthropogenically or meteorologically derived seismic signal contamination. We envision that seismic monitoring techniques will open doors to answering multiple scientific questions of the fluvial geomorphology of ice-affected streams. Apart from terrestrial streams, these results and methods can have implications for subglacial sediment transport, which are actively studied to understand changing sediment yield during climate change.
Currently, no compilation of fluvial sediment flux data exists for subarctic Fennoscandian streams. Through this study, which takes advantage of the emerging field of environmental seismology to measure bedload transport below ice and during ice break-up, we provide one of the first of such values, which is slightly higher (56.2 ± 0.7 t/km 2 /a vs. <50 or <5 t/km 2 /a) for the region than indicated by global-scale maps (Lvovich et al., 1991;Walling & Webb, 1983). Seismic techniques would allow numerous additional measurements from other catchments in the region to better constrain postglacial landscape erosion and sediment removal rates. River ice regimes are predicted to change under climate change (Prowse et al., 2011), which will likely impact catchment-scale erosion, and baseline values are required for predicting channel response and planning river management. Few studies have taken into account river morphology alterations and thus changes in bed load flux in future simulations and forecasts of northern-latitude rivers (Lotsari, Thorndycraft, & Alho, 2015). One reason for this has been the lack of verification data from all seasons, but measurements with seismometers combined with other methods will enable better verification of predictive morphodynamic models.