Modeling key processes affecting Al speciation and transport in estuaries

A R T I C L E I N F O


H I G H L I G H T S
• Al speciation is of major importance for toxicity to organisms living in estuaries. • Dynamic model estimates of estuarine transport and speciation of riverdischarged Al • High resolution model and dynamic salinity-dependent transformation of species • Model results agreed well with observed concentrations of size fractionated Al. • Key uncertainties linked to currents, model Al input and specie transformations

G R A P H I C A L A B S T R A C T A R T I C L E I N F O A B S T R A C T
Assessments of the impacts of aluminium (Al) to aquatic organisms in estuarine waters have suffered from the lack of available models that can accurately predict the presence of toxic physico-chemical forms (species) of Al at adequate spatial and temporal resolution. In the present work, transport and distribution of river-discharged Al species through changing environmental conditions in the Sandnesfjorden estuary, South-Eastern Norway, was predicted using a numerical model system at relatively high spatial (32 m × 32 m in horizontal) and temporal (1 h) resolution. New model code was implemented, including dynamic, salinity-dependent speciation and transformation processes, based on in situ measurements from several Norwegian estuaries as well as experimental data. This is the first time such elemental speciation code including LMM, colloidal, particle and sediment species is utilized in an estuary case in combination with high resolution hydrodynamics and compared to an extensive observational dataset. Good agreement was obtained between modeled and observed total and fractionated Al concentration at several stations along the fjord transect. Without including background contribution of Al from the coastal water, the model predicted too low Al concentrations (by up to approximately a factor 4) near the fjord mouth. The surface Al concentrations were also underestimated due to overestimated near-surface vertical mixing in the hydrodynamic model.
The observed correlation between salinity and total Al concentration was well reproduced by the model

Introduction
The input of trace elements and radionuclides to the marine environment is largely attributed to river transport via estuaries to the coastal zone. In estuaries, here defined as water bodies with riverine fresh water input and free connection to the open ocean, the environmental conditions are continuously changing when the fresh river water mixes with the saline coastal water. Along the estuarine transport pathways, river-discharged trace elements such as aluminium (Al) can appear in a series of physico-chemical forms (species), ranging from single ions and low molecular mass species (LMM) to larger colloids and particles (Salbu, 2009). The exact borderlines between the specie categories are difficult to distinguish as fractionation techniques are needed to separate species according to size, and thereby are operationally defined. Here, particles are defined as entities with nominal diameter greater than 0.45 lm, colloidal species are within the range from 10 kDa (few nm) to 0.45 lm and LMM species are less than 10 kDa (Salbu, 2009). In most rivers in the southern part of Norway, pH is around 6 or less, and a high fraction of Al can be present as positively charged LMM cation species (Teien et al., 2005;Schartau et al., 2016). To improve the water quality and then reduce the loss of aquatic biodiversity, for instance to avoid depletion of the wild salmon populations, most acidic rivers are limed (pH close to 7), so that Al in the river water discharging into estuaries is predominantly present as neutral species, colloids and particles. During low flow, LMM and colloidal species would be predominant, while during flooding, especially the fraction of particles will increase significantly due to erosion.
The concentration levels as well as the specie distribution and the transport properties of trace elements in estuaries are affected by a number of complex dynamic and biogeochemical processes, initiated by shifting environmental conditions such as pH, salinity and temperature, as previously observed in dynamic mixing zones (Rosseland et al., 1992;Teien et al., 2004;Teien et al., 2006b). The six most important processes considered here are: Dilution: Once the river water enters the estuary, dilution due to dispersion and mixing with coastal water masses gradually reduce the concentration levels of all river input species of trace elements being typically low in seawater. A linear negative correlation between concentration level and salinity is commonly referred to as 'conservative behavior' (Machado et al., 2016), reflecting that dilution is the predominant process influencing the concentration levels. In most cases, however, deviations from this pattern are observed due to additional processes taking place (Machado et al., 2016). Hydrolysis and polymerization: In contact with sea water, pH increases and LMM cationic species of Al will hydrolyze, polymerize and form transient colloidal positively charged polymer species being highly reactive towards available surfaces Teien et al., 2004). The colloidal fraction has been shown to contribute significantly to the estuarine transport of metals (Sañudo-Wilhelmy et al., 1996). Aggregation and sedimentation: Colloids are kept in solution due to Brownian movements, but due to increasing conductivity/salinity in estuaries, aggregation of colloids occurs. In addition, aggregates and river transported particles can settle in the estuarine sediments as the current speed decreases. Thus, the trace element concentrations of species suspended in the water column will decrease while the concentration levels in sediments will increase. Sorption, desorption and remobilization: In fresh water, a large fraction of trace elements are sorbed to surfaces of colloids and particles, i.e., associated reversibly to surfaces due to physisorption (van der Waals forces), due to electrostatic sorption (Coulomb attraction) or irreversibly due to chemisorption (covalent bonding). Upon contact with sea water, the salinity and hence the concentrations of monovalent and divalent element species will increase. Therefore, trace elements reversibly associated to colloidal and particle surfaces will be remobilized due to ion exchange processes (Teien et al., 2006b), and subsequently remobilized LMM species will be subject to hydrolysis, polymerization and formation of colloids. Remobilization of trace elements from surfaces of riverine colloids and particles as well as sediments may locally increase the total concentration and especially the LMM fraction despite high dilution effects present in estuaries (Teien et al., 2006b;Machado et al., 2016;Sanial et al., 2017). Resuspension: Particle surfaces in the sediments can act as sink for trace elements due to sorption and sedimentation. Under events with high water flow in the bottom layer and enhanced bottom stress, such as storms or high river flow rates, the particles can be mechanically resuspended and return to the water column. Formation of anions such as aluminate: Alternatively, when pH exceeds 7, LMM Al species can hydrolyze stepwise and transform to aluminate, i.e., LMM anions which is the predominant form of Al in alkaline water (pH >7.5) such as sea water (Lydersen, 1990).
Among the above-mentioned processes, the hydrolysis, polymerization, sorption/desorption and aluminate formation processes are salinity-, time-and temperature dependent, slow during winter (snow-melt episodes) and rapid during summer . Furthermore, these processes and the following speciation are also of major relevance for organisms living in the estuarine zone, as transient Al polymers formed are toxic to fish (Teien et al., 2006b). Thus, dynamic models are needed for describing the mixing zone system and processes taking place. The basis for modeling dispersion of contaminants in the ocean or the atmosphere are the three-dimensional circulation and the subsequent evolution of the density fields. This is a nonlinear initial-value problem where accurate predictions require detailed information of the initial state of the system. As was demonstrated by Lorenz (1963), due to nonlinearity of the governing equations, even very small errors in the initial conditions will grow and eventually dominate the solutions. This is popularly referred to as the 'Butterfly effect'. In weather forecasting this problem is handled by systematic acquisition of observations that are assimilated into the numerical models. In ocean prediction in general, and particularly in the present study, it is lacking both a sufficient observational and data-assimilation system for this to be possible. The evolution of predictive variables, such as salinity and temperature, is therefore only dependent on the dynamics of the numerical model itself. It is very important that the reader keep this in mind when interpreting the results of this paper.
In numerical transport estimations, where a simplified description of the real environment is required, each defined model phase (compartment) usually contains a range of species with different properties (Periáñez et al., 2018). The 'dissolved fraction' commonly includes both the LMM fraction and the colloidal fraction (Machado et al., 2016). In cases where the transport properties are different for each model compartment, such simplifications may introduce uncertainty in the model predictions (Simonsen et al., 2019). Most models applied for estuaries have focused on dilution and sedimentation, while most often processes affecting element speciation have been ignored. In addition, most experiments on trace elements or radionuclides in estuaries are based on the determination of total concentrations of filtered samples. Data sets from estuary expeditions, providing information on trace elements in different physicochemical forms, such as LMM, colloidal and particles are scarce (Lind et al., 2006). Therefore, in the present study, the overall goal was to establish a generic model system for estuaries, including LMM, colloidal and particle species and being applicable for prediction of the transport of contaminants such as radionuclides and trace element species from rivers via estuaries to the ocean, serving as input in environmental impact and risk assessments. Using the ROMS-TRACMASS numerical dispersion model system, the first objective was to implement codes describing predominant elemental species such as LMM, colloid and particles as well as key processes influencing the specie distribution in estuaries in more detail than has been done in previous literature. Secondly, the kinetics of the transformation processes were estimated based on experimental data of the behavior of aluminium species in estuaries, taking impact from factors such as salinity into account. Finally, a case study of the Sandnesfjorden estuary in south-eastern Norway was performed, where the model output was validated against field experimental data of aluminium species obtained from the river outlet and through the estuary.

Model description
The marine transport simulations in the present study were computed with the Lagrangian trajectory model TRACMASS (Döös et al., 2017), which previously has been used in a number of studies of ocean transport (e.g., Döös and Engqvist, 2007;Döös et al., 2011;Nilsson et al., 2013;Simonsen et al., 2017;2019). The three-dimensional pathways of a finite number of Lagrangian units (hereafter called trajectories) were computed 'off-line', i.e., using pre-stored input fields from an ocean circulation model. Additional horizontal diffusion was parameterized using a random walk method with constant diffusivity.
In the case study, the simulation period was from April 20, 2008 to May 31, 2008. Pre-stored three-dimensional ocean circulation fields from Sandnesfjorden were used from a hydrodynamical model simulation using ROMS (Regional Ocean Modeling System, http://myroms.org) with 32 m × 32 m horizontal resolution and one hour temporal resolution in the output. Model variables at the open boundaries were provided from a fourfold nested model system where the horizontal grid was refined from 4 km (Lien et al., 2014) to 800 m (Albretsen et al., 2011) and 160 m, all model  Fig. 3. b: The waterways upstream of the model river outlet. Dashed blue circle is the real river outlet, solid blue circle is the model river outlet. Map from www. norgeskart.no. c: Overview of the Skagerrak region, the blue square indicates the model domain. systems using ROMS, similarly to the model system applied and described by Serra-Llinares et al. (2018). Open boundary conditions were applied to transfer momentum and tracers through the open boundaries (Marchesiello et al., 2001), by providing radiation conditions on outflow and nudging on inflow. All these coastal and fjord models applied high-resolution atmospheric forcing from a 3 km simulation using the Weather Research and Forecasting model (WRF), developed by the National Center of Atmospheric Research (NCAR) (e.g., Skamarock et al., 2008). The 800 m model had a long spin-up period from January 2005, the 160 m model simulation started at March 1, 2008, and the 32 m simulation had a spin-up period from April 1, 2008 to April 20, 2008. Initial fields for the 32 m model were interpolated from the most coarse configuration.
Daily estimates of volume fluxes for six rivers (where River Storelva is the largest) in the domain were based on estimates from a distributed version of the HBV model (Beldring et al., 2003) provided by the Norwegian Water Resources and Energy Directorate. The freshwater discharges from River Storelva were supplied at Laget downstream of Lagstrømmen (the narrow strait between Naevestadfjorden and Sandnesfjorden, Fig. 1).
In the TRACMASS configuration used in the simulations in the present study, each trajectory is associated with one out of six model compartments, representing different physico-chemical forms of Al species. Aiming to represent and model the real marine behavior of metals, the model compartments and transformations were implemented as sketched in Fig. 2, where the interactions between the model compartments are implemented according to the stochastic method described by Periáñez and Elliott (2002), with transfer rates controlling the probability to shift to each of the other compartments during a time step. The trajectories can appear as two different LMM species (cations and anions) representing species with diameter less than 10 kDa, two colloidal species (humic colloids and polymers) as well as two forms of metal species reversibly bound to solid matter; seabed sediments and suspended particulate matter (SPM). The colloidal compartments represent Al species with diameter ranging from 10 kDa to 0.45 lm, either as reversibly adsorbed to the surface of suspended organic material (humic colloids) or as polymers originating from hydrolysis of the LMM species. The size fractions were selected based on operationally defined nominal sizes using fractionation techniques.
While the LMM and colloidal species flow passively with the water masses, the particle compartment represents species with diameter larger than 0.45 lm that are sufficiently large to be affected by gravity and small enough to stay suspended in the water column for significant time. In addition to the advection and turbulent mixing with the water masses, the particle-associated trajectories will experience an additional settling velocity w s towards the seabed, determined by Stoke's law (Stokes, 1851); where q s is the particle density, q w is the density of the water, g is the gravitational acceleration, d is the particle diameter and m is the water viscosity. In the model, the applied particle diameter is one single value, assumed to be representative for the average particle size. When particle-associated trajectories reach the bottom of the deepest water layer, they are transferred to the sediment compartment where they are preliminary deactivated from further advection. Due to mechanical stress from near-bottom currents, Al in the seabed sediments may also resuspend when the current speed in the deepest model layer exceeds a critical threshold value. In such case, the trajectory is reactivated as a particle-associated specie located in the deepest model layer.
The concentration levels in each of the element species are computed and displayed either as time series of the sum of trajectories present within a box around the stations within a representative time interval, or as smoothed two-dimensional fields of the surface (upper 1 m) concentration.

Sandnesfjorden case study
The present study considers estuarine transport of riverdischarged Al from River Storelva using the ROMS-TRACMASS marine numerical dispersion model system during a period with extensive Al monitoring surveys, utilizing realistic time series of river input and three dimensional hydrodynamic currents and salinity fields. The model domain covers Sandnesfjorden, an 8 km long, relatively shallow and narrow fjord located at the southeastern coast of Norway (Fig. 1). The tides are minor with maximum amplitude less than 50 cm (www.sehavniva.no). Upstream of Sandnesfjorden, River Storelva enters Songevann and Naevestadfjorden, two closed basins connected to Sandnesfjorden by Lagstrømmen, a canal which is 20-100 m wide and 3 m deep at the shallowest. Under low flow conditions, the salt wedge penetrates beyond Laget into Songevann and Naevestadfjorden basins, with salinities above 25 psu, while the surface water is relatively fresh and dominated by river water, especially at high flow (Tjomsland and Kroglund, 2010). The hydrography in Sandnesfjorden depends therefore both on the river flow and the mixing conditions. The largest depth in Sandnesfjorden is 65-70 m, while the threshold around Store Furuøy in the outer part of the fjord is around 30 m deep. Outside Store Furuøy, Sandnesfjorden is connected to the Skagerrak, a marginal sea to the North Sea, where the Norwegian Coastal Current runs predominantly southwestward along the coast. The water masses in Skagerrak are a mixture of relatively salt water from the North Sea and relatively fresh water from the Baltic Sea and adjacent coastal regions and estuaries along the Skagerrak. The dispersion simulations were performed with relevant parameters as listed in Table 1 and shown in Fig. 2b, considered to be the best possible estimates for the actual scenario.

Estimation of the transfer rates
In TRACMASS, key processes for Al speciation were computed as dynamic processes with specie interactions based on kinetic transfer rates, as sketched in Fig. 2. Since pH, ion composition and temperature in the estuary depend on the mixture between fresh water and salt water, we use salinity as the key controlling factor in this study, and the representation of the different environmental conditions from the river outlet to the open ocean was simplified by defining four discrete salinity intervals; 0-1 psu, 1-10 psu, 10-20 psu and above 20 psu. First, a set of transfer rates was suggested for each of the four salinity intervals. These suggested values were based on observed Al specie distribution from 2007 data from water samples and controlled mixing experiments at a broad range of salinities and pH values from Sandnesfjorden as well as from other estuaries such as River Lona (Teien et al., 2006b;Skalsbakken, 2009;Kroglund et al., 2011a;Kroglund et al., 2011b). The rates were estimated aiming to obtain reasonable element specie distributions after short time (< 1 day) of mixtures and after assumed equilibrium conditions after 10 days of mixing. The time evolution of concentration levels in each specie in an idealized closed system was obtained by numerical integration of the differential equations involving the interactions between the LMM, colloidal and particulate species. Interactions with seabed (sorption/desorption and sedimentation/resuspension) were neglected during estimation of the transfer rates since the sorption processes are limited to the bottom layer and the sedimentation rates depend on inhomogeneously distributed external conditions such as flow properties, settling velocity and water depth. Finally, the different rates were evaluated and adjusted against each other to obtain reasonable time scales for each process, based on the current understanding of estuarine transformation kinetics obtained from experience from numerous field studies and experimental field work (e.g. Teien et al., 2004;2006b).
These estimated transfer rates were utilized in the Al transport simulation, and predicted concentration levels were compared to observed surface concentration of the LMM, colloidal and particle Al fractions obtained at four stations in Sandnesfjorden (Hope, Skåttholmen, Hopestranda and Store Furuøy, see Fig. 1) at different times during the simulation period.

Source term description
In the model simulations in the present work, the flow through Lagstrømmen was considered to be the input of Al trajectories, released in a single grid cell at Laget in position 58.682 • N, 9.071 • E, distributed between the LMM cations, humic colloids and particle species (consistent with the observed distributions at Laget, Table 2). The water masses were initially uncontaminated. Background Al concentration was not included and no other sources than River Storelva were considered. Trajectories crossing the open boundaries of the model domain were permanently deactivated. It was assumed that the Al concentration in River Storelva depends on the flow rate similarly to that of River Lona (Teien et al., 2006b), and hence the number of trajectories released in particle compartment (N p ) increased exponentially with increasing river flow (Q); where k is a constant and the scale factor x ensures that the total number of trajectories released in the simulation is reasonable. As an approximation of the real discharges, the input concentrations of the LMM and colloidal species were assumed to be constant, i.e., the number of released numerical trajectories in these model compartments (N lmm and N c , respectively) was linearly proportional to the daily river flow, described by the following equations: The factors f p , f c and f lmm control the distribution within each compartment. The constants were set to f p = 0.2, f c = 0.5, f lmm = 0.3, x = 50 and k = 1.6. Time series of the number of input trajectories, concentration levels and specie distribution in the simulation period are shown in Fig. S.1. The magnitude of the total discharge ranged from 110 to 140lg L −1 which is comparable to the measured surface concentration at Laget (Table 2). In total, 251, 874 trajectories were released with 4-hour time intervals.

Observational data for model validation
For the purpose of validating the hydrodynamic model, measurements from 51 repeated CTD casts from Sandnesfjorden taken within the simulation period (Kroglund et al., 2011c) were used to evaluate the corresponding model output directly by calculating differences.
Surface Al concentration levels from the dispersion model predictions were validated against observational data based on surface water samples collected at stations Hope, Skåttholmen, Hopestranda and Store Furuøy (see Fig. 1 for locations), at May 10, 12 and 21, 2008, described by Kroglund et al. (2011c) and here presented in Table 2. The samples were fractioned in situ and the results are presented as concentration distribution between LMM, colloid and particle fractions after total Al concentration was determined using ICP-OES according to the method described by Teien et al. (2006b).

Kinetics of the transformation processes
To represent the natural behavior of Al in estuarine environments, our estimated rates for transformation of species (Fig. 2b, Table S.1, with corresponding time scales in Table 3) were assumed to change with salinity according to the methods described in Section 2.2.1. Estuaries are in most cases in non-steady state and in situ observations therefore reflect unstable water qualities. Unfortunately, field measurements of transfer rates are not available, but experimental data from controlled mixing experiments can be utilized, demonstrating the changes in element speciation as function of pH, salinity and time (Teien et al., 2004;2006b;Kroglund et al., 2011b,c). As the model simulations were performed during spring, the transfer rates were calibrated according to intermediate temperatures (∼10 • C). At higher temperatures, during summer and fall, we can expect faster transformations, while at lower temperatures during winter, the transfer rates will be lower (Lydersen, 1991).
The model input of Al from River Storelva was assumed to be distributed between LMM cations, humic colloids and particle fractions   (Kroglund et al., 2011b). Therefore, in the lowest salinity interval (0 psu to 1 psu), of which the purpose was to represent the distribution in the river water, we assumed that the LMM anion and polymer fractions were negligible and hence the transfer rates for the interactions including these species were set to 0. In the alkaline brackish water, mobilized LMM cations from river transported colloids and particles will polymerize rapidly to polymer-associated species when the salinity exceed 1 psu (Teien et al., 2004;2006b). Since the pH in alkaline sea water is above pH 8 (Table 2), the polymerization of LMM cation was assumed irreversible and was described by the transfer rate k 16 which was set to 0 in the lowest salinity interval, while it was high and slightly increasing with increasing salinity. LMM cations can also sorb to the available surface of solid matter present as humic colloids and SPM, controlled by the transfer rates k 12 and k 13 , respectively. These sorption rates were lower than the polymerization rates (k 16 ) and were assumed to decrease with increasing salinity due to stronger competing effects of ions in salt water and reduced available sorption surfaces of aggregated particles. Sorption of LMM cation to seabed sediments (k 14 ) was assumed to be negligible and was set to 0. Electrostatic sorption processes are reversible, and the corresponding remobilization was described by the rates k 21 (desorption from colloids), k 31 (desorption from SPM) and k 41 (desorption from seabed), respectively, which all were assumed to increase with increasing salinity. While a significant fraction of Al is assumed fixed to slowly reversible or irreversible sites at mineral particles, a large fraction of the colloidal fraction is reversibly associated to colloids in river water that easily remobilize in sea water. Hence, the rate k 21 (desorption of LMM cations from humic colloids) was relatively large and higher than the desorption rates from particles (k 31 ). The desorption rates from seabed (k 41 ) was reduced with a factor 0 compared to k 31 since the interaction surface is smaller for seabed particles than for suspended matter (Periáñez et al., 2018). Due to the combination of relatively high fraction of humic colloids in the river input and that desorption from humic colloids (k 21 ) was higher than desorption from SPM and seabed sediments (k 31 and k 41 , respectively), the input to the LMM cations was mainly controlled by the humic colloid concentration. This in turn controlled the fraction available for the rapid polymerization (k 16 ).
The polymers can dissolve as LMM anions or aggregate to SPM, which both are rapid processes. Polymerization and remobilization of LMM anions from the polymer species were described by the rates k 56 and k 65 , respectively, the latter was around 30 times higher than the first one. As these rates are sensitive to pH, they were both assumed to be unchanged with respect to salinity, as the pH was reported to increase only from 7.3 to 8.3 (Table 2).
The rate of particle formation by aggregation of humic colloids (k 23 ) is substantial, but assumed to be lower than that of remobilization (k 21 ). Due to aggregation of colloids and sedimentation of particles, lower concentration is observed in estuarine water than predicted by dilution (Hydes and Liss, 1977), and the aggregation was assumed to increase with increasing salinity due to suppression of the electrical double-layers of particles (Lebovka, 2012). Similarly, the rate describing aggregation of polymers (k 63 ) was assumed to increase with increasing salinity. The effect of mechanical stress on the particles resulting in decomposition (weathering) to smaller species (k 32 and k 36 ) was considered negligible and hence these rates were set to 0.
The interactions between LMM and the solid matter (SPM or seabed sediments) have in some ocean transport studies successfully been simulated as two-step functions, with one compartment where the trace element is reversibly bound (adsorbed) to the surface and one with slowly reversible/irreversible bindings (e.g. Periáñez, 2003Periáñez, , 2004. These processes are however expected to play a minor role in short-term and near-shore simulations (Periáñez et al., 2018) and were therefore ignored here.

Model assessments
To validate the model hydrography and predicted Al concentration in the present study, data obtained from a number of CTD records and water samples have been used. With high spatial and temporal resolution, the current model configuration provides high detail levels in the results with obvious advantages, especially when local small-scale processes with rapid fluctuations are investigated. Validation of model results with in situ measurements may, however, be challenging. Due to the fundamental underlying differences between the concentration levels achieved from numerical models and from in situ observations, some discrepancies are expected when comparing model results with observations (Sandvik et al., 2016). While models compute average values in a finite volume and time period, observation samples are snap-shot values attributed to a certain point. In regions with non-uniform distribution, this representation error associated to the samples may therefore be considerable (Janjić et al., 2018). Although the resolution in our model is considered to be high in most contexts, small-scale natural processes at sub-grid scales, that cannot be properly resolved by the model, may be significant here as well, especially close to the river outlets, and need therefore to be parameterized. Other factors that may affect the model performance are smoothing of the coastline and topography in the model, response time in the model and negligence of background concentration of Al in inflowing coastal water. However, our study is aiming to describe the Al speciation and transport through the whole fjord, at relatively high spatial and temporal resolution, utilizing a suggested set of estimated transfer rates. Based on transformation kinetics (Section 3.1), the transfer rates were chosen aiming to achieve reasonable distributions at relatively short as well as long time scales utilizing local measurements where applicable. Hence, to resolve the details of the speciation and transformation Polymer → P a r t i c l e -1 2 4 . 6 3 . 5 processes accurately, high resolution was required. Although the very rapid transformation processes (e.g., minutes) were omitted, this makes the model able to provide concentration levels in the complex speciation setup through the estuary from the river outlet to the open ocean.

Validation of hydrography in the hydrodynamic model
The comparison between the measured CTD-profiles and corresponding values from the hydrodynamic model revealed that the model in general reproduced the observed salinity and temperature relatively well. However, the model predicted too high salinities (7 psu in average) in the surface layer (Fig. S.3). The salinity error decreased linearly with depth and was close to zero at about 13 m below surface, indicating that the offshore salinities were realistic. The temperatures in the model were around 2 • C too low in the surface layer and around 1.5 • C too high at depths below 11 m.
The model bias in near-surface salinities was not evenly distributed along the fjord. Comparison between the measured surface salinities and the corresponding model predictions show large variability and a relatively small bias at the stations in the inner part of the fjord, while the model overestimated the salinities closer to the fjord mouth ( Fig. S.4). The first indicates that the applied freshwater discharges were realistic, although there may be deviations according to timing of the run-off pulses. The latter indicates that the modeled near-surface salinity exaggerated the alongfjord gradient, and this may be attributed to too weak vertical stratification and deficiencies in the vertical mixing in the ocean model.

Strong wind event May 12, 2008
Strong winds will have a huge effect on the vertical positioning of the pychnoclines along the shore, either by elevating the pychnoclines during offshore winds (upwelling) or lowering during onshore winds (downwelling). The distribution of the water masses may then be altered for hours or days. Fig. 3 shows predicted along-fjord transects 1 of modeled salinity for the upper 20 m before, during and after a strong wind-driven overturning of the water masses occurred in the morning of May 12, 2008. The onshore winds reached gale force from the north-east 2 . In the model, this event induced inflow of coastal water masses into Sandnesfjorden and subsequent deep mechanical mixing with downwelling of the near-shore water masses, then resulting in an increased nearsurface salinity in the outer part of the fjord. Before the event (May 11, 2008 at 21 UTC, Fig. 3a), the isohalines were mostly organized horizontally with a thin low salinity surface layer. The freshwater input from River Storelva is clearly seen as low salinity water in the inner part of the fjord. At the end of the strong wind period (May 12, 2008 at 09 UTC, Fig. 3b), the outer part of the fjord became more vertically well-mixed, with weaker stratification and higher surface salinity. At the same time, the model shows that due to the strong wind-driven surface currents pushing coastal water into the fjord, the discharged river water was blocked in the inner part of the fjord, and thereby reducing the salinity through the whole water column near the river outlet. The isohalines became more vertically aligned than before the wind increased. In the afternoon, when the wind had ceased (May 12, 2008 at 15 UTC, Fig. 3c), the accumulated river water flushed out of the fjord, resulting in a low salinity layer near the surface as the water masses returned towards the original state.
These rapid exchanges of water masses seen in the model suggest that the exact location, timing and sampling depth will have large impact on the measured values. With strong gradients, the representation uncertainty in the observations is considerable, and similarly, small displacements of physical features will affect the model results significantly. Model-observation comparison under such circumstances has to be done carefully, since small spatial or temporal displacement of a local phenomena can lead to large deviations in the comparisons.

Trajectory age at the stations
In order to demonstrate how fast the trace elements were transported away from the release point at Laget to different locations in Sandnesfjorden, histograms of the modeled age of trajectories present in the surface water at the four stations are shown in Fig. 4. Most of the trajectories arrived at Hope less than six hours after they were released, with the peak value between 2 and 3 h. The age increases with distance from the source, and at Store Furuøy, near the fjord mouth, the peak was between 24 and 36 h. Hence, as the time resolution in the hydrodynamic fields as well as the phase shift was 1 h, there were few possibilities for the model trajectories to obtain the correct phase distribution already at Hope, even with high transfer rates. A better agreement with observed distribution of Al species could therefore be expected at the stations at the greatest distance from the source.

Total concentrations
Horizontal distribution of modeled total Al concentration (sum of all Al species) is shown in Fig. 5, for a date with high river flow early in the simulation period (Fig. 5a) and for dates later in the period, when the observation samples were taken (Fig. 5b-d). Due to dilution of river water, the estimated concentration levels generally decreased with distance from the river outlet, with occasional accumulation of Al in near-shore bays and coves along both shores. While neglecting background level of Al in the coastal water, the modeled concentration levels decreased towards zero in the far-field regions. Due to generally decreasing Al discharges through the time period, the Al inventory in the fjord also decreased with time. A clear effect of the strong wind event is seen at May 12 (Fig. 5c), with remarkably higher surface Al concentration in the inner part of the fjord and reduced Al concentration levels in the outer part compared to the other dates.
Included in the figure are also surface observations of total Al concentration from Hope, Hopestranda and Store Furuøy at May 10, May 12 and May 21, as well as from Skåttholmen at May 12 (Table 2). In general, the model results agreed well with measured data, although the Al concentrations were underestimated in the outer part of the fjord at all times, primary attributed to the negligence of background contamination in the coastal water. This is confirmed by the time series in Fig. 6, showing time series of estimated surface total concentration and salinity from the model simulation as well as surface samples of observed total Al concentration and salinity (Table 2) and the surface values (upper 1 m mean) of salinity from the CTD casts taken at the stations. As previously shown (Section 3.2.1), the modeled surface salinity was systematically overestimated at all stations (7 psu in average), mostly near the fjord mouth and in the beginning of the period. Later in the simulation period, the agreement with observed surface salinity was pretty good (<5 psu) at all stations.
At least in the outer part of the fjord, the dilution of river water, causing increased salinity and decreased surface concentration as the riverine water propagate through the fjord, was overestimated by the model. This points to a too weak stratification and too strong vertical mixing of the water masses in the hydrodynamic model. This is a well-known feature related to coastal modeling, especially under stable stratified conditions (Zilitinkevich et al., 2007). Numerical models tend to smooth out such strong gradients due to discretized model layers and imperfect parameterizations of turbulent flow features. Such large gradients caused by small-scale processes can therefore not always be expected to be properly resolved by the models. It appears that in this particular case, the water masses were strongly stratified with a thin surface layer consisting of almost pure river-water. This layer was only slightly diluted by the coastal water, and hence high Al concentration and low salinity was observed even in the outer part of the fjord. This was supported by the observed Al concentrations which hardly decreased from Hopestranda to Store Furuøy. The water masses in the model were, however, more effectively diluted due to a combination of too strong upwelling of uncontaminated deep water and too strong vertical turbulent mixing, and the concentration levels decreased with around a factor 5 between Hopestranda and Store Furuøy (Fig. 5). This resulted in overestimated salinity (>10 psu) and underestimated Al concentration (up to factor 4), especially in the outer part of the fjord. In our discharge scenario, only Al releases from River Storelva were considered. Generally, other potential sources that in reality may contribute to the observed metal contamination are releases from nearby rivers, diffuse surface run-off such as agriculture and roads, atmospheric deposition and discharges through the groundwater (Machado et al., 2016). Particularly in our case, output from nearby rivers mixing Al into the coastal current upstream (north) of Sandnesfjorden were expected to contribute to some extent. Assuming constant Al concentration in the coastal water, the contribution of background concentration should be linearly correlated with the salinity. Hence, to account for Al concentration in surrounding Skagerrak water, a background term was added to the model predicted concentration (C), with the relation: whereĈ was the background corrected concentration. C 0 and S 0 were set according to data taken at approximately 20 m depth in Flødevigen at May 22, 2008 (Table 2), considered to be representative for the background levels in surrounding coastal water. Time series of estimated concentration levels inclusive this background term are shown with the black lines in Fig. 6. As intended, the Al concentrations increased mainly in the outer part of the fjord, while it was almost unchanged in the inner part. Although this background correction could not completely account for the deviation between the model results and observations at Store Furuøy, the inclusion of the background level of Al improved the model output.
Looking at the other available background data (Table 2), the variability in observed Al concentration was disproportionately high, despite small variations in salinity. Hence, according to our limited dataset, the deep water cannot be assumed to keep constant background concentration, and other additional, and presently unknown, important factors were causing the variability. Even though there are reasons to believe that a substantial fraction of the Al contamination in the outer part of Sandnesfjorden came from other sources than River Storelva, applying a salinity-dependent background term as in Eq. (5) will introduce large uncertainties associated with both algorithms and input data. Therefore, the results in the rest of this paper will be presented exclusive the background term, although the reader should be aware that the results should principally be corrected for this term.
Single occasional episodes such as the strong wind event at May 12 (discussed in Section 3.2.2) also affected the model results heavily at all stations, with increased mixing of deep-water in the outer part of the fjord. In the inner part of the fjord, a notable difference between model results and observations is seen. While the model predicted accumulation of river water with decreasing surface salinity and increasing Al concentration, this was not evident in the observations. Actually, at all stations, even at Hope and Laget (Table 2), the observed surface salinity increased considerably (>5 psu) between May 10 and May 12, while the Al concentration decreased with more than 30%. This indicate that in this particular situation, the entrainment of deep-water to the surface was stronger in reality than predicted by the model, even in the inner part of the fjord, close to Laget. Neither was the predicted accumulation of river water seen in the observations. In addition, processes taking part upstream of Laget, outside the model domain and hence not captured by the model, might also have contributed to disagreements between model and observations.
Numerous other factors could also have contributed to discrepancies between model and observations. As discussed above, a slightly temporal or spatial displacement of the dynamical structures may result in large errors when comparing to a few in situ samples. From the model results, we have identified large spatial gradients as well as rapid exchanges of the water masses. Therefore, the representation errors may be considerable, and we cannot assume all the estuarine dynamics and variabilities to be captured by the observation samples. Due to strong vertical stratification and near-surface gradients, the measurements were affected by the exact depth of the sample equipment, where small perturbations in the sampling depth, caused by waves and motions in the vessel, might have affected the results. The measured salinity can thus be expected to be slightly underestimated.In addition to the factors regarding the flow field from the hydrodynamic model, model errors can also be attributed to uncertainties in parameterizations and assumptions in the configuration of the transport model such as uncertainties in the discharges.

Conservative behavior
Due to dilution in estuaries, river-discharged contaminants with low sea water concentration are expected to show so-called ' conservative behavior' with a negative linear correlation between concentration and salinity (Machado et al., 2016). Such a relationship was observed in Sandnesfjorden (Fig. 7a), where the field experimental Al concentration and salinity data taken during our simulation period were strongly negatively correlated (Pearson r = −0.97, Table 4). Despite these correlated observational data, when extracting model estimated time series from fixed positions, the temporal variability of Al concentration and salinity were only slightly correlated. The correlation was low at all stations (Table 4), highest at the stations near the river outlet (Hope and Skåttholmen), and close to zero at Hopestranda and Store Furuøy. As expected, the model variability was higher than the observed, due to a higher number of ' sample pairs' from the model (n = 1004) than in the observations (n = 10), covering a broader range of conditions. The salinity bias and negligence of background concentration were both expected to affect the results as discussed in Section 3.3, but due to the large uncertainties involved, these corrections were not included in the model results shown in the present subsection. In contrast to the low correlation at each station, the correlation in total for all stations was relatively high (Pearson r = −0.80, Table 4). This indicate that the model predicted a general trend of conservative behavior when considering Sandnesfjorden as a whole. However, some deviations from this conservative behavior occurred, especially in the results from Hopestranda (Fig. 7a). Although the Al concentration in general decreased and the salinity increased with increasing distance from the river outlet, occasionally, both the Al concentration and the salinity increased from Hope to Hopestranda. As seen in the maps in Fig. 5, Al was accumulated with high concentration gradients close to the shore near the location of Hopestranda. Only a gentle displacement of these fronts in the model could have large impacts on the results. Hence, the high Al concentrations predicted at Hopestranda might have been attributed to such small gradient disagreement between model and reality. Deviations from the conservative behavior are also expected to be caused by uncertainties in the model representation of current flow and mixing properties, variabilities in salinity in the coastal water and uncertainties in the Al concentration and freshwater flow rate in the river. In addition, transformation of Al species will complicate the transport due to sorption/desorption processes as well as aggregation, precipitation and resuspension of colloidal and particle species.
With wind forcing, freshwater input, tides and baroclinic pressure gradients being the most important contributors for controlling the surface currents in the fjord, we expected the surface flow pattern to be a major contributor for the deviation from conservative behavior. Therefore, to identify times when the model predicted conservative behavior, in contrast to situations in which the results were non-conservative, the model results were divided into two groups depending on the upper layer volume flux 3 through a cross section between Hope and Skåttholmen, marked with a blue line in Fig. 1. The first group consisted of situations at low flux conditions, while the high flux conditions were in the second group, shown as time series in Fig. 8a. The figure also shows the bias in the surface salinity (comparison of corresponding model and CTD data) computed in each of the four stations.
While the magnitudes of both the flux and the surface salinity bias were highest in the first days of May, there was a general decreasing trend through the remaining time period for both variables. The surface salinity bias and the flux were slightly correlated (r = 0.5) with lowest bias under low flow conditions, thereafter increasing bias with increasing volume flux (Fig. S.5). The first group, consisting of times with low volume flux (below average, 33 m 3 s −1 ), is indicated with green colors in Fig. 8a, while times with high volume flux (above average, 33 m 3 s −1 ) were in the second group, indicated with red colors.
These two groups were recognized as different regimes in the water masses, affected by river flow rates and wind forcing, as shown by the Hovmöller plot of the salinity field from Skåttholmen site in Fig. 8b. Situations in the low flux regime were characterized by relatively low river flow rates (Fig. S.1) and prevailing northeasterly winds (not shown), accumulating freshwater in the inner part of the fjord. In the period between April 29 and May 2, the thickness of the surface freshwater layer exceeded 2 m depth (Fig. 8b). In the period after May 12, there were moderate winds from shifting directions and weaker stratification of the water column. The high flux events appeared in time periods with weak winds from changing directions (but mostly from southwest) and high river flow rates. The salinity increased rapidly between the surface and the well-mixed saline water below the halocline which was located closer to the surface than it was during low flow. Thus, the strong flux may be responsible for upwelling of saline deep water to relatively shallow depths.
Although variations in the flux cannot fully explain the occasional non-conservative behavior in the model, the two different flux regimes (Fig. 7b-c) were clearly distinguishable. The low flux events (Fig. 7b) appeared to be in good qualitative agreement with the observed pattern, with a clear conservative trend, i.e., high negative correlation between salinity and Al concentration (Table 4). In contrast, the high flux events (Fig. 7c) were characterized by increasing Al concentration between Hope and Hopestranda and lower correlation (non-conservative behavior).
This behavior may also partially be explained by assumptions related to the model configuration, such as resolution. Although horizontal resolution of 32 m is considered to be relatively high in most oceanographic contexts, the results were most likely affected by factors such as un-resolved sub-grid scale dynamics, transformation of species and by the model displacement of the river outlet. The non-conservative behavior could also have been caused by local near-shore features at Hopestranda, and by imperfect model descriptions of vertical mixing and upwelling processes. Combining these factors, the model might have overestimated the surface current velocities, especially in the inner part of the fjord. Too strong currents near the discharge point would have overestimated the flushing of Al away from the river outlets. Further out, where the current speed became lower, the Al contaminants would accumulate, giving rise to increasing Al concentration between Hope and Hopestranda. However, due to a limited number of field experimental results, (samples were only collected from three days during the simulation period), it is not clear if the occasions of non-conservative behavior shown by the model under high flow events were a part of the natural variability or if they were caused by artificial or numerical mis-predictions introduced by the model system. But what is clear is that these occasional events where the model estimated increasing Al concentration from Hope to Hopestranda appeared almost exclusively during high flow events. Hence, according to the available observations, the model performed best during low flow.

Distribution of Al species
Contour plots of modeled surface concentrations in each of the elemental species are shown in Fig. 9, for a date early in the simulation period with high river flow (May 2, 2008, Fig. 9a-c) and a date later in the simulation period with lower river flow (May 21, 2008, Fig. 9d-f). At the May 21 plot, the model results are compared with observed surface concentration of Al species. The pattern of the LMM species ( Fig. 9a and d) was quite homogeneously distributed through the fjord with relatively high concentration levels (40-50 lg L −1 ) reaching almost out to the mouth of the fjord. In the inner part of the fjord, the LMM fraction consisted mainly of cations, while the outer part was dominated by anions which were dissolved from polymers (not distinguished in figure). The colloidal (Fig. 9b and e) and particle species (Fig. 9c and f), in contrast, had the highest surface concentration levels in the inner part of the fjord, which decreased towards the fjord opening and thus had stronger gradients along a fjord transect than the LMM fraction. Since the particles settled towards the bottom, they were effectively removed from the surface water layer, especially during low flow (Fig. 9f). Under high flow, the surface particle concentration was relatively high also in the outer part of the fjord (Fig. 9c). Scatter plot of estimated Al concentration and salinity for each specie in Fig. 10 was in good agreement with our intended behavior of the transformation processes (Section 3.1). With river discharges distributed between LMM cation, humic colloids and particles, the highest concentrations of these species appeared at low salinity at Hope, near the river outlet. Affected by dilution, the concentration of these species decreased with increasing salinity. However, due to rapid remobilization of Al reversibly associated to humic colloids and subsequent polymerization of LMM cations, the concentration of these species decayed faster than the linear dilution, dependent on the time of the reactions. The polymers were formed from the LMM cations at high transfer rates (k 16 ), and were assumed to be unstable with short life time. Hence, the pattern of the polymer species was similar to that of the LMM cation species, except at the very lowest salinities. The concentration of LMM cations and polymers of Al decreased rapidly with increasing salinity, and both were highest near the river outlet. Predicting these species is essential since they are assumed to be the most toxic to biota (Teien et al., 2006a). As the remobilization of polymers was fast (k 65 ), while the polymerization of LMM anions was slower (k 56 ), the life time for the LMM anions became relatively long. Therefore, as a consequence of higher mobilization rate than dilution close to the river outlet, the concentration of LMM anion species increased with increasing salinity from Hope to Hopestranda, reaching a maximum at around 15 psu. By correcting for the salinity bias (Section 3.2.1), the modeled LMM anion concentration would have reached maximum at around 5-10 psu, in agreement with Kroglund et al. (2011c). Further out in the fjord, from Hopestranda to Store Furuøy, the LMM anion concentration decreased due to dilution. The prediction of LMM anions is also of interest as these species are known to be reactive towards fish gills (Teien et al., 2011).
Time series of estimated and observed Al specie distribution at the stations are shown in Fig. 11. The observed colloidal fraction was close to zero and was overestimated by the model at all stations, especially at Hope and Skåttholmen where the model predicted ∼30% colloids. At the same time, the LMM fraction was underestimated by the model, also this mostly at the innermost stations. At Hope, the model predicted ∼40% LMM species, while the observed fraction was ∼75% (Fig. 11). This may be explained by the negligence of background Al in coastal water (mainly present as LMM species, Section 3.3), the assumptions regarding specie distribution in the discharges (Section 2.2.2) and the relative short distance between the river outlet and the first stations (Section 3.2.3). With only a few hours travel time after release, and hence only a few possibilities for transformation of species, the trajectories did not have sufficient time to obtain the expected specie distribution. This issue Fig. 10. Scatter plot of estimated surface Al concentration in each specie vs salinity. Each circle represent one model 'sample pair', obtained at 4 h time intervals, from the four stations marked in Fig. 1a. was however improved in the model results at Hopestranda and Store Furuøy, the stations located most far away from the source as well as in the late time period at all stations.
Comparison with observed fractions shows that the particle fraction was relatively well predicted by the model at all stations. The response of the river flow is clearly seen in the time series for the predicted particle fraction at Hopestranda, as it increased at high flow and decreased under low flow conditions. Higher concentration of Al particles further out in the fjord under high flow conditions was also seen in the maps in Fig. 9. An opposite pattern can be seen for the LMM fraction. Hence, at locations around Hopestranda, it can be assumed that the particle species mainly originated directly from river discharges, requiring high river flow rates and strong surface currents out of the fjord to reach that station. In addition, there are large uncertainties related to the role of interactions with the sediments. Sedimented particle-bound Al can be remobilized to become LMM cation species or be subject to resuspension to become particle species. Both processes will increase the transport time compared to the fraction that had not interacted with the seabed.
The combination of the results presented above indicate that to capture all the physico-chemical processes near the river outlet, the applied transfer rates might have been too slow. Alternatively, it may indicate that the travel time from river discharge to the trajectories reached the first stations was too short in the model simulations. In addition, the Al discharges were in reality most likely affected by transformation processes during the transit through the basins upstream of Lagstrømmen (Songevann and Naevestadfjorden, Fig. 1), and by sub-grid scale processes in the mixing zone which were not Fig. 11. Time series of estimated and observed Al specie concentration and Al specie fractions in surface water (upper 1 m mean) at the four stations marked in Fig. 1a. explicitly resolved by our model. Hence, increasing the model resolution may be expected to improve the model skill in the inner part of the fjord. However, these model-observations comparisons demonstrate that the model in general was capable to reproduce the observed concentration levels of Al species further out in the fjord. Although some weaknesses related to the above-mentioned issues in the current model setup limit the prediction skill somewhat with respect to toxicity, especially directly outside the river outlets, the model can be used to predict distribution of Al species at locations and times of interest within the model domain. Regarding the toxicity for living organisms, the speciation of the trace elements and radionuclides is essential (Teien et al., 2006b). Hence, utilization of such model transport simulations including dynamic transformation of species are highly needed and relevant for evaluation of the elemental specie distribution and other conditions in different parts of an estuary. By utilizing the link between elemental speciation, bioavailability, uptake in organisms and toxicity, such a model can be used for risk assessments.

Conclusions
Numerical transport models are useful and applicable tools for predictions of total concentration and metal speciation in the marine environment over large areas and over continuous time periods. The results provided here are unique as this is the first time transport and dynamic transformation kinetics involving colloidal species of a trace element in a fjord system has been predicted over a considerable time period, utilizing high resolution three-dimensional hydrodynamic model fields (32 m horizontal resolution). The present study is also unique as it compares predicted distribution of Al species to real observations at several locations and times during the simulation period.
In general, results from the simulations validate relatively well with observed data of Al fractions taken in Sandnesfjorden in 2008. However, the quality of the results are affected by a number of more or less uncertain assumptions and input data involved in the model configuration. The observed correlation between salinity and total Al concentration (conservative behavior) appeared to be best predicted by the model under low flux conditions. Although the Al concentrations in the study area were densely monitored, with fractionation data taken at a several times and at a number of stations through the Storelva-Sandnesfjorden estuary, large uncertainties are still related to the description of the transformation processes. In the present case, imperfect hydrography predictions might have affected the transport estimates, and discrepancies between modeled and observed total Al concentration in the fjord mouth indicate significant contribution of Al from the coastal water, probably discharged from rivers through adjacent fjords and coastal areas. To reduce the uncertainties associated to the background levels in future monitoring fieldwork, we suggest that the Al concentration in coastal waters upstream of the key area are properly assessed, ideally by taking samples at different depths in surrounding waters during the whole monitoring period. The model predicted too weak vertical hydrography gradients, and hence, too strong mixing of fresh river water with the more saline coastal water was assumed to affect the model performance as the large near-surface gradients appeared to be smoothed out by the model. Other factors that were assumed to affect the results are imperfect implementation of the river discharges and parameterization of sub-grid scale processes. Hence, we suggest even higher temporal and spatial resolution for further improvement of the model system, as well as a better description of the processes involved during transit through the river basins upstream of Lagstrømmen.
However, despite these issues, as the general trends of Al concentration and speciation were well reproduced, the results obtained by present work demonstrates the potential of high resolution models with dynamic speciation to predict exposure of contaminants such as metals and radionuclides for risk assessments in the marine environment. This is especially relevant for models that are able to provide predictions of the bioavailable species at high detail level in the full water column, not only in the surface. As the tool we have developed and applied in the present study is generic, it can easily be extended for application of other metals and radionuclides.