The case study site constitutes the contaminated fiberbank outside of Väja pulp and paper factory. The fiberbank is situated in the inner part of the Ångermanälven river estuary, along the northern Bothnian Sea coast, Sweden (Figure 2). The estuary is shaped as a fjord-like basin that formed in a prolonged tectonic rift (Heinemo, 1999). The mean river discharge is ∼500 m³ s^–1 at the river mouth. Between the years 1999 and 2015, a maximum of ∼2,300 m³ s^–1 and a minimum of ∼90 m³ s^–1 were recorded. The water in the estuary mostly contains 70–80% freshwater, but at discharges above 830 m³ s^–1, the inflow of seawater is blocked and the proportion of freshwater temporarily becomes higher (Heinemo, 1999).

During the peak discharge of oxygen consuming organic fiber material from the industry in the post-war years of the 1950s, a clear expansion of anoxic seabed in the Ångermanälven river estuary was identified (Jerkeman and Norrström, 2018; Valeur, 2020). Many industries have closed down since then (Heinemo, 1999; Apler et al., 2020), and a recovery has been observed in some of the areas. However, several contaminated fiberbanks have been identified along the river estuary (Apler et al., 2014, 2019), and the mere presence of a fiberbank means that it has affected the natural sediment substrate and habitat at that location.

The studied contaminated sediment site is located adjacent to a sulfate pulp and paper factory (Figure 3) that has been in operation since the 1910s. Up to the 1980’s, oxygen consuming organic fiber material had accumulated on the seabed close to the factory’s water outlet. Although the discharge of suspended solids (including fibers) has essentially stopped, a large fiberbank, divided into two parts, remains, together covering an area of ∼70,000 m² (Apler et al., 2019). The fiberbank is anoxic and consists mostly of cellulose fibers with a thickness exceeding 6 m in places. The fiberbank is resting upon natural minerogenic sediments some 20 m below the water surface and the surface of the deposit is located at ∼15 m water depth (Apler et al., 2019). The natural sediments successively slope about 10 degrees seaward down to a water depth of ca 50 m some 200 m out from the shoreline. The seabed next to the bank contains fiber-rich sediments that have been estimated to cover 800,000 m², i.e., >10 times larger than the fiberbank itself. Several scars from submarine landslides are visible in 3D models of the seabed created from bathymetrical data (Figures 3, 4). The factory quay today receives about 40 vessels a year (personal communication with staff from Väja factory) for the handling of goods.

All field measurements are presented in Apler et al.(2014, 2019, 2020), Dahlberg et al. (2020), Dahlberg et al. (2021), and Löfroth et al. (2021).

Geological and geotechnical surveys (2015 and 2017) show that the fiberbank is underlain by a thin layer of bark on top of a clayey silt, followed deeper down by laminated black-gray clayey silts that transition toward the base into glaciomarine varved clay (Löfroth et al., 2021). An interbedded coarse-grained layer was detected about 1–2 m down into the minerogenic sediment. The geotechnical investigation indicates that this layer has higher permeability and higher porewater pressure than the surrounding sediments (Löfroth et al., 2021).

Field samples were taken during a field study in Väja in August and September 2015. Subsequent laboratory studies of the samples showed that the fibers easily swirl during physical disturbance and that the fibers tend to tangle up and form small clumps (Lamparski, 2016). This observation shows that the resuspended fibers are thus of sufficiently low density to stay in suspension, disperse by currents and not settle until they reach calmer areas outside the fiberbank. Further, laboratory studies on particles obtained from two fiberbanks in the study area also indicate low settling velocities and indicate that the density of the pulp fibers is similar to the density of water (Löfroth et al., 2021). Figures 5, 6 show two photos of the fiberbank and fiberbank material at Väja.

Sediment-to-water fluxes of PCBs were measured in situ using benthic flux chambers (BFCs) equipped with semipermeable membrane devices (SPMDs) as passive samplers (see Eek et al. (2010) for detailed description of the BFCs used in this study). Seven indicator PCB congeners (CB-28, −52, −101, −118, −138, −153, −180) were measured and their sum concentrations (Σ₇PCB) were used to illustrate the suggested methodology. The following matrices were analyzed with respect to levels of Σ₇PCB: sediment (fiberbank, ng g^–1 dw); porewater (pg L^–1); diffusive flux (ng m^–2 day^–1) from SPMD samplers placed on the seabed; and dissolved fraction of the PCBs in sediment pore water (%). In addition, sediment (fiberbank) organic carbon content (% TOC), was determined. Detailed information on the chemical analyses is given in Dahlberg et al. (2020), Dahlberg et al. (2021).

River flow and wind data were obtained from SMHI. These data were used to theoretically calculate possible current conditions that may affect fiber particle mobility.

Dispersion pathways considered in this study were chemical diffusion and mechanical advection caused by current induced resuspension (propeller wash, wind-waves, river flow velocity) and submarine landslides (Figures 7A–E). The risk can be explained as a chain of events. In our study, the risk chain contains the following three events:

If the probability of event 1 is equal to 1.0 (100%), the probability of exceeding a relevant contaminant threshold can be estimated by combining event 2 and 3. Here, event 2 concerns dispersion through five separate dispersion pathways and exceedance is estimated for each of them. Event 2 is expressed as P_X in Eq. 1 and event 3 is expressed as P[F_X > I_T] in Eq. 1, further explained below.

Assessing the likelihood of contaminant dispersion thus consists of two parts: (1) assessing the likelihood that each of the five dispersion pathways in Figures 7A–E will occur within a given time; and (2) assessing the likelihood that one or several dispersion pathways will result in contaminant dispersal above a certain impact threshold considered likely to have negative impact on the environment. A general formula for assessing the probability of impact follows Göransson et al. (2014), Göransson et al. (2018) and is

where P is the probability, P_I,X is the probability of impact (I) for a particular dispersion event (X), P_X is the probability that this particular event will occur within a given time (here set to 1 year), F_X is the flux of contaminants from such an event if it occurs, and I_T is the impact threshold. For each of dispersion event, a method to calculate P_x and F_x is determined. For the calculation of F_x, parameters are designated a probability distribution to account for uncertainties in the parameter values. Fluxes are presented in the unit of mg m^–2 year^–1.

In our study, we measured diffusion (Figure 7A) in the field using BFCs (Dahlberg et al., 2021). For other pathways (Figures 7B–E), we must rely on mathematical calculations based on sediment and contaminant characteristics measured in the field or using field samples, and empirical relationships found in the literature. For the specific dispersion pathway, X is here replaced by D for diffusion, S for ship (propeller) current induced resuspension, Q for river flow current induced resuspension, WC for wind-wave current induced resuspension, and L for submarine landslide-induced resuspension.

The method presented in this paper has been applied on contaminated, organic (cellulose) rich sediments of anthropogenic origin but can also be applied to contaminated sediments in general. The method aims to support risk assessment by estimating the risk of dispersion, considering, but not limited to, dispersion through chemical diffusion and advection induced by particle resuspension. We have chosen to use a probabilistic approach, as a complement to deterministic risk assessment, to consider systematic and statistical uncertainties with the aim to assess the likelihood of pollution spread and its likelihood for environmental impact. However, one difficulty with such an approach is the quantification of the overall uncertainty due to a forward propagation of uncertainties, as well as model and parameter uncertainties (e.g., Gouldby et al., 2010). On the other hand, a probabilistic method may reveal uncertainties that otherwise can be hidden if only a deterministic analysis is applied (Arabi et al., 2007; Qu et al., 2016). As such, a probabilistic method may be more rigorous, for example by displaying percentiles and confidence intervals that are more precise (Qu et al., 2016) and considering different environmental conditions, as shown by Shojaeezadeh et al. (2020). Although the confidence intervals here can be considered large, they give an understanding of possible best- and worst-case scenario for the different continuous, periodic or rare occurring dispersion pathways based on the available knowledge. For example, by comparing the Monte Carlo simulated fluxes and loads (Table 3) with calculated fluxes and loads for each of the dispersion pathways by using fixed mean and maximum values (except for Manning’s M for which the minimum value should be used to yield higher shear stress) (Supplementary Table 4), we found a discrepancy for river flow current, wind wave current, and submarine landslide induced dispersion. In the case of limited data, calculations based on individual values may therefor give misleading results.

Using a probabilistic approach helps to understand where the uncertainties lie and identify parameters for which uncertainty needs to be reduced to better estimate the risk. Such information can form a basis for the design of supplementary field studies and laboratory analyzes.

The method presented in this study provides information on whether or not dispersion occurs, but does not predict where the dispersed contaminants and particles will end up. For such an assessment, a numerical model is needed. Toxicological screening and modeling could also be included to assess the impact on the ecosystem being exposed. Thus, the accuracy to which the risk can be estimated or assessed with the suggested method depends on how well the threshold level for environmental impact can be described. Nevertheless, a probabilistic method can be used to confirm the validity of a deterministic risk assessment, if such exists (Shojaeezadeh et al., 2020).

In our study, the computational model was set up in a spread sheet tool, supported by Monte Carlo simulation, which makes it easy to update the (mathematical) analytical solutions when new knowledge become available, if governing equations need to be changed, and if additional dispersion pathways need to be added. Also, Monte Carlos simulations assume variables to be independent. However, some of the uncertain independent variables may not be truly independent. For example, diffusion depends on concentration gradient, wavelength depends on wind speed and fetch, near bed flow velocity depends on discharge and bed geometry. When more data are available, one could choose assumptions to correlate in a correlation matrix with correlation coefficients as input to the sensitivity analysis.

Based on the results, it can be concluded that Σ₇PCB are dispersed and that fiberbank particles are resuspended from the studied fiberbank. The results clearly show that the rare occurrence of extreme events will have a relatively large impact on the amount dispersed, which indicates that the fiberbank cannot be considered stable with respect to contaminant and particle dispersion. The results further indicate that the fiberbank will likely continue to negatively affect the surrounding ecosystem if nothing is done to prevent dispersion. This outcome is probably also true for many other fiberbank deposits, and for contaminated sediments located in shallow waters (e.g., Butcher and Garvey, 2004; Davis et al., 2007; Howell and Rifai, 2015). Even archived contaminants may become exposed to dispersion because of sediment scouring during extreme weather-events like hurricanes (Dellapenna et al., 2020). Our study also includes the potential for contaminant dispersion from submarine landslides, and the results indicate that a submarine landslide will disperse most contaminants per unit area (i.e., flux), followed by advection through ship induced resuspension, river flow resuspension, wind-wave current resuspension, and least by chemical diffusion.

We applied a threshold level for Σ₇PCB sediment-to-water flux to estimate the probability for environmental impact. Such an approach was also applied by Shojaeezadeh et al. (2020). However, in our study, this threshold is related to measured Σ₇PCB fluxes at a reference site because no Swedish environmental quality standard exists for PCBs in water. The threshold level thus simply reflects what the situation is in a less contaminated part of the area, but it does not give information on potential consequences to the environment or the recipient, which would have been the better threshold level. The Norwegian Miljødirektoratet has set a guideline value for a yearly average of Σ₇PCB in coastal waters to 2.4E-06 μg L^–1 (Miljødirektoratet, 2016). Sediment-to-water fluxes can be recalculated to an approximate water concentration if the representative water volume and the water residence time is known. For the Bollsta Bay, the water volume is about 0.11 km³ (SMHI, 2021a) and the residence time is roughly 1 year, depending on water depth (Heinemo, 1999). In our case study, we chose a threshold level of 0.00013 mg m^–2 year^–1, which roughly would correspond to 4.9E-06 μg L^–1. Thus, the suggested threshold level is about two times higher than the Norwegian guideline, but within the same order of magnitude, which means that the selected threshold value is not unreasonable. The water volume above Väja fiberbank, not the entire Bollsta Bay, is roughly 1 million m³ and with the same residence time, the measured yearly fluxes of chemical diffusion from the fiberbank (0.0011 mg m^–2 year^–1) would roughly correspond to a concentration of 7.7E-05 μg L^–1, which is about 30 times higher than the Norwegian guideline value. Calculated fluxes from the other dispersion pathways are even higher than that. To conclude, this means that the estimated risk for our case study is an underestimation rather than an overestimation, even considering the lower levels of the 90% confidence intervals. It also confirms the importance of taking into account advective dispersion pathways in risk assessments (e.g., Erickson et al., 2005; Davis et al., 2007; Reible, 2014; Howell and Rifai, 2015; Dellapenna et al., 2020).

The uncertainty in input parameter values can be reduced with more precise laboratory and field measurements. We used benthic flux chambers (presented in Dahlberg et al., 2021) to measure chemical diffusion, however, we cannot be sure that we only measure diffusion with the BFCs because the advection can also be influenced by e.g., bioturbating benthic organisms. We did not have the possibility to measure all the parameters needed for the calculations, partly because existing field and laboratory equipment are not designed for the atypical material that the fiberbanks consist of. Thus, several parameter values had to be estimated based on literature and the knowledge gained so far about this anthropogenic sediment. There are some parameters that were more challenging to assess. For example, the high organic cellulose fiber content and their characteristics made it difficult to assess parameters like critical shear stresses, grain size distribution and bed surface roughness, which are parameters commonly included in formulas for sediment transport. It was also challenging to conduct proper geotechnical measurements of the fiberbank, for example to collect undisturbed core samples of the soft unconsolidated fibrous sediment and the underlying strata (Löfroth et al., 2021). The varves (thin silt and clay layers) in the underlying soft clay made trimming of specimens for laboratory tests challenging (Löfroth et al., 2021). Flow and current velocity profiles would have provided better input to calculate shear stresses. The optimal would have been to have access to measurement data showing current directions and current velocities, calibrated against prevailing flows, winds and ocean currents, and then to establish a hydrodynamic model across different scenarios.

It has also been shown that resuspension of fiberbank materials and fiber-rich sediments may be induced by only gentle disturbance at the site (Dahlberg et al., 2021), which should be considered in future work. Gas production (methane gas from anaerobic decomposition of organic matter) can also be high in fiberbanks (Lehoux et al., 2021), and gas ebullition may perhaps enhance contaminant dispersal in fiberbanks through increased sediment-pore water exchange (Klein, 2006) and through resuspension of contaminated particles (Yuan et al., 2007; Viana et al., 2012). Gas formation can cause parts of the fiberbank to be boyant and float away. Furthermore, the geological setting at a site may provide conditions for advection through upwelling submarine groundwater. In our case study, the upwelling of submarine groundwater was not assumed an issue because the permeable thin silt layer is not in contact with the contaminated sediment layers. The observed pockmarks can indicate submarine groundwater seepage but are more likely caused by gas ebullition.

In the middle and northern parts of Sweden, the effect of climate change includes increased precipitation and milder winters, which in turn will cause more precipitation to fall as rain instead of snow. A climate analysis for the county shows that the run-off to the rivers in the area increases as the climate changes (SMHI, 2021b). Thus, the flow in the Ångermanälven river is likely to increase, which may affect the flow-induced resuspension of fiber particles. However, the Ångermanälven river discharge is strongly regulated by several hydropower stations, and an ongoing landslide risk investigation of the Ångermanälven river valley conducted by SGI has shown that the flow control obscures the effect of climate change on river flow. Nevertheless, increased rainfall will affect pore water pressure and hence landslide probability. It is not currently known whether climate change will cause stronger and/or more frequent storms in the area, but if these occur, wind-wave currents will pose a greater problem than today for the dispersion of contaminants and fiber particles. Other studies have also shown that sediment-to-water dispersion of PCBs increases during hurricanes (Howell and Rifai, 2015; Dellapenna et al., 2020).

The case study site is located in a part of Sweden that currently experiences isostatic rebound of the earth’s crust of 8 mm year^–1 in response to the last glaciation more than 10,000 years ago, thus resulting in a gradual seaward displacement of the shoreline. This means that the water depth over the fiberbank is slowly decreasing, eventually exposing the fiberbank more to currents from both vessels, discharge and waves. However, this situation may change over time due to sea level rise, and it is possible that the rate of sea level rise may exceed the isostatic rebound rate in about 100 years (personal communication, L. Johansson, SMHI), but by then, the sediments may already have been redistributed over the seabed. In addition, isostatic rebound and other tectonic movements of the earth’s crust cause minor earthquakes, but so far the highest measured magnitude in the area is 3.1 on the Richter scale¹ and it has been considered that these magnitudes are not enough to trigger a landslide in the area (Löfroth et al., 2021).

The scientific knowledge on the characteristics of the cellulose fibers and wood chips deposits in the environment is in its infancy. The fiberbank deposits will be difficult and costly to remediate due mainly to high organic content, high water content, very low density, and the large volumes (Jersak et al., 2016a,b). It may also be difficult to determine who is legally responsible for the costs and implementation of the action in accordance with the Polluter Pays Principle (PPP). Because not all sites can be remediated at once, we need to understand the impact from these deposits on the ecosystem in order to prioritize sites for remediation. We need to understand more about the different dispersion pathways from the fiberbanks to the surrounding recipient and about their different significance for the spread of these pollutants. This can lead to an understanding of the fiberbanks’ impact on the ecosystem and in turn how we address fiberbanks to minimize their environmental impact. The approach presented here takes uncertainties into account and can be used as an initial assessment of whether dispersion can be a problem, which dispersion pathway probably contributes most, and in the dispersions order of magnitude.