{'fraction_non_alphanumeric': 0.03436310943753327, 'fraction_numerical': 0.01934033870236382, 'mean_word_length': 5.189627821842587, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling Reactive Solute Transport in Permafrost‐Affected Groundwater Systems

Understanding the interactions between ground freeze‐thaw, groundwater flow, and solute transport is imperative for evaluating the fate of contaminants in permafrost regions. However, predicting solute migration in permafrost‐affected groundwater systems is challenging due to the inherent interactions and coupling between subsurface mass and energy transport processes. To this end, we developed a numerical model that considers coupled groundwater flow, subsurface heat transfer, and solute transport, including water‐ice phase change, solute‐dependent porewater freezing, and temperature‐dependent solute reaction rates. As an illustrative example, we present simulations to investigate the potential for contamination from a municipal wastewater lagoon in the Canadian sub‐arctic. Two‐dimensional groundwater models assuming varying permafrost conditions were developed to evaluate possible contaminant migration scenarios associated with groundwater flow from the lagoon to a river, including the transport of conservative, degradable, and sorbing solutes. Model results reveal important transport mechanisms controlling the behavior of aqueous contaminants in permafrost landscapes as well as the hydrogeologic factors affecting reactive transport in cold regions. Seasonal freeze‐thaw episodically restricts connectivity of transport pathways, which attenuates both transport and reaction rates. However, elevated solute concentrations can depress the freezing temperature of porewater and produce thaw‐induced solute transport. Both thermally driven and solute‐enhanced thaw can decrease ice content in permafrost, which can have significant implications for solute migration. Therefore, thermo‐hydrogeologic process controlling reactive transport in cold‐region groundwater systems must be considered when quantitatively assessing the impact of changing environmental conditions on contaminant hydrogeology.

surface warming, the possibility exists for dissolved solutes from these facilities to leach into groundwater and groundwater-fed waterways, potentially posing a risk to ecosystem health, and into drinking water supply relying on surface water (Burkow & Kallenborn, 2000;Hayward et al., 2018;Krumhansl et al., 2015;Macdonald et al., 2000). Considerable effort is required to identify, delineate, and remediate subsurface contamination. The use of groundwater flow and transport models to support management efforts is common practice; however, to address the potential for groundwater contamination in cold regions there is a need to adapt solute transport descriptions in mathematical models to account for the presence of pore ice and freezing temperatures (Panday & Corapcioglu, 1994;Vonk et al., 2019). Present solute transport models generally do not consider permafrost dynamics, while cold-regions groundwater models including solute transport routines are rare (Grenier et al., 2018;Lamontagne-Hallé et al., 2020). Previous research on subsurface contaminant transport in cold regions has predominantly focused on petroleum spills and active-layer transport of other non-aqueous phase liquids (Carlson & Barnes, 2011;Dyke, 2001;Panday & Corapcioglu, 1994). Few studies have focused on the behavior of aqueous solutes within active-layer and permafrost groundwater systems on relatively small spatial scales (Bosson et al., 2013;Frampton & Destouni, 2015;Yi et al., 2021). Consequently, there is a need to improve understanding of how changing permafrost conditions may impact groundwater flow and reactive solute transport dynamics in cold regions.
Although saturated ice-rich permafrost has a very strongly reduced permeability (up to seven orders of magnitude) compared to when liquid water occupies pore space, perennial taliks and active layer groundwater flow during seasonally unfrozen periods provide pathways for advective and dispersive contaminant migration above and below permafrost (Bosson et al., 2013;Johansson et al., 2015) as well as diffusion into permafrost containing unfrozen pore water (Barnes & Chuvilin, 2009;Biggar et al., 1998). Seasonal freeze-thaw and associated weathering processes can also enhance active-layer hydraulic conductivity and increase solute advection and dispersion (Ding et al., 2019;Koch et al., 2013;Mohammed et al., 2021). Additionally, the degradation and transformation of reactive solutes are driven by temperature-dependent biogeochemical activity and freezing conditions that mediate the fate of reactive solutes (Dusek et al., 2019;McCarter et al., 2017;Yurova et al., 2008). Improved conceptual and mathematical descriptions of coupled groundwater flow and reactive solute transport processes specific to cold regions are thus necessary to assess potential contamination in permafrost groundwater systems (McKenzie et al., 2021). Such model development is also needed for investigating changing groundwater-facilitated transport dynamics brought on by climate change and permafrost thaw (Bense et al., 2009;Frampton et al., 2013;Ge et al., 2011;Kurylyk et al., 2016;Walvoord et al., 2012), such as the transport of nutrients and heavy metals released from thawing permafrost soils (Neilson et al., 2018;Selroos et al., 2019;Vonk et al., 2019).
In this study, we developed a numerical cryohydrogeological model that considers groundwater flow, porewater freeze-thaw, and reactive solute transport using physics-based descriptions of subsurface water, heat, and chemical transport, including constitutive relationships for water-ice phase change, sorption, and temperature-dependent solute degradation. The new model will be useful for studying how permafrost conditions impact groundwater flow and reactive solute transport for different contaminants across a range of geologic, hydrologic, and climatic conditions. As an illustrative example, we apply the model here in the first instance to investigate potential groundwater contamination associated with a wastewater lagoon facility in the Canadian subarctic. These simulations identified important transport mechanisms controlling the behavior of conservative, sorbing, and degrading aqueous contaminants in permafrost regions. Results are interpreted to support site investigations to assess the possible impacts of the wastewater lagoon on nearby groundwater and surface water quality.

Model Development
The generic finite-element modeling environment FlexPDE (PDE Solutions Inc., 2020, http://www. pdesolutions.com) was used to solve the partial differential equations and constitutive relationships describing coupled groundwater flow, reactive solute transport, and heat transfer with freeze-thaw. Cryohydrogeological models built within FlexPDE have been previously validated against benchmarked numerical solutions of groundwater flow and permafrost thaw (Grenier et al., 2018) and applied to permafrost groundwater systems without solute transport (Bense et al., 2009(Bense et al., , 2012. This section describes the water, heat and reactive transport equations used in this study. Assuming saturated subsurface conditions, water flow is described by (Bense et al., 2009): is acceleration due to gravity, μ [kg m −1 s −1 ] is dynamic viscosity of water, k [m 2 ] is the intrinsic permeability of porous medium, k rw [-] is a relative permeability exponential function that reduces permeability as pore ice is formed (McKenzie et al., 2007), and S s [m −1 ] is the porous medium specific storage. Energy transport and subsurface temperatures (T [°C]) are modeled using the conduction-advection equation considering the effects of latent heat associated with the melting/freezing of pore ice: where κ e [W m −1 K −1 ] and C a [J m −3 K −1 ] are respectively the effective thermal conductivity and volumetric heat capacity of the sediment/fluid/ ice mixture, θ w [m 3 m −3 ] is the liquid water content, and  q [m s −1 ] is the Darcy flux (k rw kh), i.e., groundwater flow in Equation 1 drives advective heat transport. C w [J m −3 K −1 ] and L i [J m −3 ] are the volumetric water heat capacity and volumetric latent heat of fusion, respectively. The reduction in effective liquid water saturation, S w [-] as temperatures decrease and ice forms is described using an empirical exponential function approximation of the soil freezing curve (SFC) (Lunardini, 1988;McKenzie et al., 2007): where S wres [-] is the residual liquid saturation and W (K) is an empirical parameter related to the temperature transition window over which liquid water decreases from saturation to S wres as ice forms (Lunardini, 1988).
Experimental evidence has shown that certain soluble salts (e.g., Chloride (Cl − )) decrease the freezing temperature of porewater by increasing the osmotic potential of the aqueous solution in the porous medium (Cary & Mayland, 1972;Ginot et al., 2020). The net effect of this phenomenon on the SFC is to increase the amount of liquid water present at sub-zero temperatures and increase the temperature window of the water-ice transition period or the "slushy zone" of the SFC (Ginot et al., 2020). In this study, solutes are not expected to significantly affect porewater freezing beyond this effect (as maximum measured Cl − concentrations in the lagoon were approximately 100 mg L −1 , and seawater Cl − concentration is around 19,400 mg L −1 for reference). Thus, here we employed a simple modification to the SFC to consider the effects of concentration-dependent freezing as follows. First, the SFC for groundwater with a Cl − concentration of 100 mg L −1 was approximated using the salinity-dependent SFC formulation by Zhou et al. (2020). The curve fitting parameter W in Equation 3 was then adjusted to approximately match the change in SFC (Figure 1). Thus for Cl − solute transport, W is linearly increased from 0.5 to 1.8 over the solute concentration range of 0-100 mg L −1 (Figure 1). This formulation enables us to represent the effect of elevated Cl − concentrations on the liquid water present and thus the permeability for temperatures from 0 to −3.5°C (Biggar et al., 1998;Panday & Corapcioglu, 1994). Expressions for all constitutive relations and parameters values used in Equations 1-3 are listed in Table 1.
In addition to the mediating effects of pore ice on advective-dispersive transport and concentration-dependent porewater freezing, the reactivity of solutes in groundwater involves temperature-dependent degradation processes such as microbially driven mineralization and biodegradation, and physio-chemical processes such as adsorption/desorption (Fan et al., 2010). To make the model numerically tractable when coupling solute transport to both water flow and heat transfer, simplifying assumptions had to be made. (a) The effects of solute concentration on fluid density were not considered (see above discussion of lagoon Cl − concentrations). (b) An equilibrium approach was used to describe solute sorption (Dusek et al., 2019).
(c) First-order kinetics were assumed for solute degradation (Yurova et al., 2008). Solute reaction processes include sorption and degradation as mass sources or sinks in addition to advective-dispersive transport: where C [mg L −1 ] is the solute concentration, D [m 2 s −1 ] is the hydrodynamic dispersion tensor defining solute spreading by mechanical dispersion and molecular diffusion whose terms have been modified for the presence of pore ice (Table 2), v i [m s −1 ] is the linear fluid velocity ( accounting for the presence of pore ice (Table 2), with the groundwater velocity from the flow equation (Equation 1) driving the advective solute flux. Terms ρ b [kg m −3 ] and n [m 3 m −3 ] are the bulk density and porosity of the porous medium, and S C [g kg −1 ] is the adsorbed concentration of solutes onto the aquifer matrix. Solute reactions are coupled to heat transport through the degradation term, λ a [s −1 ], which is the temperature-dependent first-order coefficient for solute degradation. The temperature dependency and associated coupling to heat transport in the degradation term is described using a modified Van't Hoff equation describing the exponential increase in the rate of biological activity with temperature (Yurova et al., 2008 Volumetric latent heat of fusion for water, L i J m −3 Ice density × mass-based latent heat = 3.03 × 10 8 (Bense et al., 2009) Residual saturation, S Wres -0.05 W K 0.5 ("clean" water) to 1.8 ("dirty" water) Liquid water content, θ w m 3 m −3 θ w = n S w Hydraulic conductivity, K w m s −1 K w = k rw k ρ w g/μ Impedance factor, Ω -50 where λ a ref is the reference degradation rate for the reference temperature, T ref [°C], and Q 10 is the factor for increase in metabolic rate per 10 °C increase in temperature. Sorption, S C (Table 2) is modeled assuming a linear equilibrium sorption approach (Dusek et al., 2019). The implementation of Equation 4 in FlexPDE was benchmarked against the classic Ogata and Banks (1961) analytical solution, and Equations 1-3 have been extensively benchmarked in Grenier et al. (2018). Equations 1, 2 and 4 are solved simultaneously so that solute, heat, and groundwater transport are coupled using constitutive relations in Tables 1 and 2. Expressions for all relations and values for all parameters used in Equations 4 and 5 are listed in Table 2.

Illustrative Model Application
Our new model was applied to investigate the potential for groundwater contamination associated with a wastewater lagoon facility in the Northwest Territories, Canada ( Figure 2). For this site, models of groundwater flow, subsurface heat transfer, and reactive solute transport were developed and parameterized based on measured data including local climate, vegetation cover, soil and sediment properties, topography, groundwater depths, and seasonal air temperature variations. However, the lack of permafrost and active layer thickness data at this stage, along with the high probability of thawing permafrost in the area, precluded the ability to do a thorough data-model comparison or calibration (e.g., Albers et al., 2020). Thus, the specific objective of this modeling exercise was to evaluate the potential for subsurface solute migration associated with groundwater flow from the lagoon to reach the river for a range of conditions. Simulations were run for different scenarios of initial permafrost distributions (i.e., surface temperature conditions) to generate a range of possible contamination transport pathways associated with the water flow into the aquifer from MOHAMMED ET AL.  Bear, 1972) Vertical dispersion component, D yy (Bear, 1972) Longitudinal dispersivity, α L m 10 Transverse dispersivity, α T m 1 Factor for increase in metabolic rate per 10 °C increase in temperature, Q 10 -1.7 (Yurova et al., 2008) Reference temperature, T ref°C 5.0 Aquifer matrix bulk density, ρ b kg m −3 1,250 Initial concentration for conservative solute mg L −1 100 Initial concentration for total N (degrading) mg L −1 80 Initial concentration for total P (sorbing) mg L −1 10 Table 2 Parameter Values and Expressions Used in the Solute Transport Equation the lagoon and its impact on the development of a contaminant plume in the underlying groundwater and its potential to exfiltrate into nearby surface water.

Field Site
The hydrogeological system underlying a wastewater lagoon (Figure 2a), is located near a small rural community in northern Canada, with a population of approximately 110 people. The exact location of the community cannot be specified due to confidentiality reasons. The lagoon ( Figure 2b) is located approximately 900 m northwest of a river in the region. The region has an average annual precipitation of 336 mm, with average winter air temperatures of −22°C in January and +16°C in July (Government of Canada, 2019), with an average annual air temperature of −2°C. However, permafrost in the region, where present, is relatively warm and generally closer to −1°C (Smith et al., 2009). The site is located in the subarctic zone of discontinuous permafrost (Heginbottom & Radburn, 1992;Smith et al., 2009), and is expected to be underlain by permafrost to some degree. As is the regional trend, the subsurface at the site is likely warming, triggering permafrost thaw. Currently, the active-layer groundwater flow season has a duration of approximately 4-5 months, from May to September/October, but is possibly year-long where permafrost has already degraded substantially and a supra-permafrost aquifer has developed.
The lagoon is approximately 30 m in diameter and 2 m deep. The lagoon water surface is consistently 1-2 m below the ground surface and is connected to the local water table (Figure 2b). The lagoon is surrounded by a mixture of trees and forested areas and marshy/wetland areas, indicating shallow/local groundwater discharge areas within 500 of the lagoon. Municipal wastewater is delivered to the lagoon daily. Site investigations in 2016 show that the topography of the site slopes toward the southeast, draining toward the river ( Figure 2). Lateral hydraulic gradients measured via monitoring wells (installed to a depth of 6 m below ground surface) in the vicinity of the lagoon ( Figure 2a) were on average approximately 0.01 m/m, and indicated groundwater flows east toward the river. For the purposes of this modeling investigation, the shortest flow path to the river (Figure 2a), was assumed to be the dominant groundwater flow direction to simulate a worst-case contamination scenario for the river. Borehole logs and soil/sediment surveys generally classified the geology at the site as glacial till consisting of sandy gravel or sandy clay underlain by clay-rich till. Slug tests in monitoring wells indicated that hydraulic conductivities in the upper sandy gravel and sandy clay units where permafrost is absent ranged from 10 −4 to 10 −6 m s −1 respectively. Field investigations also identified several water quality parameters (e.g., ammonia-nitrogen, chloride) that were above background concentrations in the monitoring wells (

Permafrost and Active-Layer Properties
The cyro-hydrogeological model developed here can incorporate spatiotemporal distributions of hydrothermal properties in multiple dimensions, but is currently limited to two-dimensional mass and energy transport in the saturated subsurface for numerical tractability. Thus, estimates of ground and near-surface temperatures are needed to provide upper thermal boundary conditions (Evans et al., 2020;Kurylyk et al., 2016;Langford et al., 2020). Permafrost and active-layer depths are governed by subsurface heat transfer, which is driven mainly by thermal conduction between the ground surface and underlying subsurface (Hayashi et al., 2007;Sjöberg et al., 2016). Ground surface temperatures are determined by the surface energy balance of heat fluxes due to radiation, sensible heat, latent heat (i.e., evaporation, condensation), and ground conduction (Hinzman et al., 1998).
To make use of available meteorological data and account for these factors, the surface energy balance of the site was initially simulated using a one-dimensional soil-vegetation-atmosphere-transfer model which accounts for ground freeze-thaw conditions. The SHAW model (Flerchinger & Saxton, 1989;Mohammed et al., 2017) was used to estimate the expected ground temperatures and active-layer and permafrost-table depths. The initial estimates of ground temperatures and active-layer depths from the SHAW simulations were used to inform the range of surface thermal boundary conditions employed in the groundwater model. These were then combined with available site geology, topography and groundwater information to develop the domain and boundary conditions for the groundwater model. Information regarding the SHAW simulations can be found in the Supporting Information.

Model Domain, Initial Conditions, and Boundary Conditions
The groundwater model domain employed represents a two-dimensional topography-driven groundwater flow system underlain by permafrost. Topography at the lagoon site slopes gently toward the river to the southeast ( Figure 3). A 2 m-resolution digital elevation map of the site was used to extract the surface topography along a 2-D transect that represented the minimum distance between the lagoon and the river ( Figure 2a). As the water table is relatively shallow at the lagoon facility, it was assumed to follow surface topography. Hydraulic head along the top of the domain was fixed to represent a water table boundary condition ( Figure 3a). Consequently, recharge occurred at the ridge of the hillslope, and regional discharge was focused in the topographic low adjacent to the river (Figure 3d). The base and vertical sides of the model were "no-flow" hydraulic boundaries to represent catchment divides (Figure 3b). For heat transport, a uniform temperature of 2°C was assigned to the base of the domain (i.e., below the permafrost base) rather than a geothermal heat flux (Ge et al., 2011) to form permafrost from the top-down, while the vertical sides were closed for heat exchange (Figure 3b). The upper specified temperature boundary condition varied according to the simulation scenario (discussed below). Initial conditions shown for groundwater flow and permafrost conditions ( Figure 3) were developed by running the model to a dynamic equilibrium condition (simulated for 1,000 years) with an average annual temperature of −1°C (or −2°C for one scenario discussed below) applied to the top of the domain to generate permafrost. The generation and impact of taliks are discussed below. As FlexPDE employs time-step control and dynamic adaptive mesh refinement, element size and vertical discretization ranged from 0.01 to 10 m during simulations, with finer mesh discretization needed in the active layer and permafrost.
MOHAMMED ET AL.  To simulate contaminant transport associated with lagoon effluent, a specified concentration boundary condition for the contaminant of interest was applied along a portion of the water table under the lagoon (Figure 3a). For other regions along the top boundary, a boundary condition was specified such that any recharging water had a concentration of 0 mg L −1 . The concentration of water exiting the domain, that is, groundwater discharge, was determined by the simulated advective solute flux leaving the model, representing the possible contaminant loading to the river. The base and sides of the model were closed to mass fluxes (Figure 3). The model initially had no contaminant mass present in the system. Three contaminants were simulated: (a) a conservative solute; (b) total nitrogen, a degrading solute; and (c) total phosphorus, a sorbing solute.

Model Scenarios
Five simulation scenarios are presented as illustrative examples of possible permafrost, groundwater flow, and solute transport conditions (Table 3). For the first simulation scenario, referred to as S1, the top domain temperature boundary condition was specified such that the seasonally maximum active-layer depths were ∼2-2.5 m in the groundwater model. The temperature along the top boundary was prescribed using a sinusoidal seasonal temperature cycle with a mean of −1°C and amplitude of 8°C to force a seasonal active layer and drive episodic shallow active-layer groundwater flow (Figures 3b and 3d). Simulations were run from initial conditions for 50 years without (S1a) and with (S1b) the modified concentration dependent SFC (Figure 1). Scenario S2 (S2a and S2b) were similar to S1 but had permafrost generated from an average surface temperature of −2°C. A seasonal signal with a mean of −2°C and amplitude of 8°C was then applied to generate colder and, therefore, thicker permafrost and thinner active layer conditions. Simulations S3 and S4 were run with continuous surface warming (thawing permafrost) and the presence of taliks with no modification to the SFC. S3 has a talik present under the lagoon, generated by increasing the initial average surface temperature over the lagoon from −1°C to 5°C. S4 has taliks under both the lagoon and river, with the lagoon talik generated in the same manner as S3 and the river talik generated by applying an average surface temperature of 2°C over the river. The difference in lagoon talik configuration between S3 (closed talik in Figure 3c) and S4 (open talik in Figure 3d) is due to the presence of the river talik which increases the hydraulic connection between the lagoon and river. The implications of this connection are discussed later in the results section. These simulations were subjected to a constant linear surface warming signal of 0.05°C/year for 100 years (consistent with warming scenario RCP 8.5 (Collins et al., 2013)) from their initial condition.
Scenario S5 represented a highly idealized model domain with flow restricted to a thin, shallow aquifer and no seasonal freezing or permafrost present to maximize the rate of contaminant transport and examine the "worst-case" scenario for contaminant loading to the river. Scenarios S1 to S5 were designed to evaluate the possible range of permafrost conditions and the effect on contaminant loading to the river (Table 3). Simulation scenarios S1 and S2 were run for 50 years, S3 and S4 for 100 years, and S5 for 500 years. Specified concentrations applied to the lagoon boundary conditions (Table 2) were based on water samples taken from the lagoon. The parameters for contaminant degradation and sorption were estimated based on available literature for similar processes and/or systems, and the temperature dependency in terms described in Section 2 (Dusek et al., 2019;Yurova et al., 2008). Values for all parameters are listed in Tables 1 and 2.

Results and Discussion
The model results presented in the following sections are discussed and interpreted to elucidate the complex interplay between coupled water, heat, and solute transport dynamics in permafrost-affected groundwater systems. Figure 4 illustrates the simulated hydraulic head and hydraulic conductivity distributions for scenarios S1-S4. In scenarios S1 and S2, during the active-layer flow period, the flow regime was shallow and flow rates near the lagoon (Figure 4a) were approximately 10 −7 m s −1 . Reduced hydraulic conductivity of the permafrost (Figures 4b and 4d) restricted flow and thus the potential transport of contaminants to the shallow active layer during unfrozen periods. Also, the surface topography generated shallow local flow reversals away from the regional gradient to the river (Figures 4a and 4c), which combined with the shallow flow system, seasonal freezing in the active layer, and stable permafrost at depth, limited transport from the lagoon directly to the river. Figure 5 shows the temperature distributions at the beginning and end of the simulations for all scenarios. For scenarios S1 and S2 (Figures 5a-5d), temperature conditions were stable and did not change during the simulation period.

Effects of Permafrost Conditions on Groundwater Dynamics
Results from scenarios S3 and S4 show that the presence of taliks had a significant effect on groundwater flow dynamics (Figure 4), and thus, taliks are also expected to have a strong control on the spatial distribution of solutes from the lagoon. For S3 with a talik present only under the lagoon, essentially all flow was confined to the thawing active layer in the supra-permafrost flow system (Figure 4e), as the talik's hydraulic connectivity was limited vertically ( Figure 4f); however, there was increased downward flow in the talik itself which increased temperatures below the lagoon (Figures 5e and 5f). In contrast, in scenario S4, with taliks beneath the lagoon and river, the downward flow enhanced both the hydraulic (Figures 4g and 4h) and thermal (Figures 5g and 5h) connections to the river, which vertically opened the lagoon talik compared to scenario S3 (Figures 4e and 4f) and facilitated sub-permafrost flow in addition to the supra-permafrost system.
The range of groundwater flow conditions in scenarios S1-S4 are significantly different, representing: (a) a shallow episodic active-layer flow system in S1 and S2; (b) a supra-permafrost flow system in S3; and (c) a combined supra-permafrost and sub-permafrost flow system in S4. These differing flow, thermal, and freeze-thaw regimes would thus be expected to exert a strong control on the maximum extent and distribution of solute mass present in the system. MOHAMMED ET AL.  . Left panel shows simulated hydraulic head distributions, and right panels show simulated hydraulic conductivity distributions accounting for the presence of pore ice for: (a and b) Scenarios S1 and S2 (active layer closeup during thawed period, note scale difference from panels below); (c and d) Scenarios S1 and S2 at t = 50 years; (e and f) Scenario S3 at t = 100 years; and (g and h) scenario S4 at t = 100 years.

Contaminant Transport and Migration in Stable Permafrost
Contaminant plumes and permafrost conditions at the end of the 50-years simulation period of S1 (stable climate and seasonal active layer at an average surface temperature of −1°C) are shown in Figure 6, for plumes with the unmodified SFC (S1a, panels a-d), and for the concentration-dependent SFC (S1b, panels e-h). For scenario S1a, the conservative solute plume (Figure 6b), defined as the extent of the 10 mg/L isoline, extended to approximately 5 m below the domain surface; however, the bulk of the mass was within the upper 2.5 m, coinciding with the maximum seasonal thaw-depth of the active layer. The plume migrated MOHAMMED ET AL.

10.1029/2020WR028771
10 of 20   laterally approximately 50 m toward to the river, and about 30 m in the opposite direction due to local flow against the regional gradient that is toward the river, highlighting that regional groundwater movement in the region around the lagoon was limited by seasonal freezing in the active layer and relatively shallow flowpaths, as observed in earlier studies for other sites and settings (Lamontagne-Hallé et al., 2018;Sjöberg et al., 2016). The reactive nitrogen plume (Figure 6c), defined as the extent of the 1 mg/L isoline, was far less persistent due to degradation combined with temporary immobilization from seasonal freezing. In this case, the degradation rate rivaled the advective flux of the contaminant driven by groundwater flow, resulting in minimal horizontal migration (∼10 m) of the reactive contaminant in these simulations. The phosphorus plume (Figure 6d), defined as the extent of the 1 mg/L isoline, was also less persistent than the conservative solute, as most mass was adsorbed during active-layer transport. Both nitrogen and phosphorus plumes range in depths from 2-3 m below the base of the lagoon, which was near the base of the active layer.
Simulations for scenario S1b (Figures 6e-6h), where Cl − concentrations affect the SFC (Figure 1) show that elevated salt concentrations can potentially have an important effect on contaminant migration, even at levels far below that of seawater and neglecting variable density effects on buoyancy-driven flow components. All plumes showed enhanced migration relative to their corresponding simulations in S1a ( Figure 6). In particular, the conservative salt and phosphorous plumes showed noticeable increases in mass transport due to solute-induced thaw (Figures 6g and 6h). For example, the largest lateral plume extent (chloride in Figure 6f) was increased by 27% for the scenario with a solute-dependent SFC compared to the scenario without (chloride in Figure 6b). This solute-induced thaw locally decreased pore ice around the lagoon which increased the effective permeability of the permafrost (Figure 6e) and allowed for further solute migration compared to S1a. The nitrogen plume, however, did not show an increase in mass compared to S1a (Figure 6g), as the degradation rate was still higher than the advective fluxes. These results show that when solutes allow more liquid water to be present at freezing temperatures, transport may be enhanced due to the increased permeability of permafrost as it becomes "slushy" near the thawing point (Ginot et al., 2020).  Figure 7 shows results for scenarios S2a and S2b, which had colder permafrost due to the lower average surface temperature (−2 °C). The maximum seasonal active layer thickness (Figure 7a) in these scenarios was approximately 2 m, which was only 0.5 m less than for S1 (Figure 6a). However, in scenario S2a, the chloride plume's horizontal extent (Figure 7b) was reduced by ∼150% compared to S1, but increased vertically by 40% (∼2 m) (Figure 7b vs. Figure 6b). Less horizontal transport in this scenario increased the vertical concentration gradients around the lagoon causing more downward diffusion compared to S1a. The nitrogen plume (Figure 7c) did not show a similar increase, and was only slighter larger in extent than for S1a due to the degradation rate overpowering this increase in diffusion (Figure 6c). The phosphorus plume mass ( Figure 7d) also only increased modestly with depth compared to S1a (Figure 6d), due to increased sorption.
For scenario S2b, with the solute-dependent SFC, the increase in vertical transport of chloride in S2 compared to S1, translated into a similar increase in solute-induced thaw (Figure 7e), and corresponding increase in chloride concentrations with depth ( Figure 7f). The degrading nitrogen plume again did not change much compared to S2a or S1. However, the increase in solute-enhanced thaw caused the phosphorus plume extent (Figure 7h) to considerably increase with depth compared to S2a and both scenarios in S1, despite colder permafrost conditions. Figure 8 shows the difference in ice-saturation during the winter period when surface temperatures were at −1°C (Figures 8a and 8b) and −2°C (Figures 8c and 8d). These subtle differences in S1b and S2b in Figure 8 show that the influence of solute-enhanced thaw is less important for colder permafrost conditions, that is, it becomes less noticeable as temperatures decrease below the temperature-transition window of the SFC, but this can still have important implications for the migration of solutes (e.g., chloride and phosphorus in scenarios S1 and S2). From a modeling perspective, the simplified implementation of the concentration-dependent SFC limits the description of W in Equation 3 as an effective parameter that would need to be adjusted depending on the range of salinities present at the site being modeled. More physically based relationships between salinity and soil freezing (Zhou et al., 2020) would circumvent this. Another important numerical implication of these results is that if salt concentrations do affect the SFC, then multiple solutes need to be simultaneously simulated to produce the effect of increased liquid water content and consequent enhanced permeability on the migration of other solutes Figure 7. Results for S2a and b with stable permafrost at a surface temperature of −2°C, that is, similar to S1 but under colder conditions. Panels (a-d) show results of scenario S2a at simulated time = 50 years showing (a) icesaturation distribution at the maximum seasonal active layer thaw-depth, and plumes for (b) chloride, (c) nitrogen, and (d) phosphorus. Panels (e-h) show results for scenario S2b (modified SFC) with (e) ice-saturation distribution at the maximum seasonal active layer thaw-depth, and plumes for (f) chloride, (g) nitrogen, and (h) phosphorus. present in the system. S1b and S2b were ran with all three species being simulated in parallel, which significantly increased run times.

Contaminant Transport in Thawing Permafrost
Scenarios S3a and S4a were run from a stable initial condition (surface temperature = −1°C) and forward modeled subject to a constant surface warming of 0.05°C/year so that thaw continuously occurred (Figures 5g and 5h). At the end of 100 years in S3a, with a closed talik under the lagoon, the conservative plume was pooled at the base of the thawed layer with little diffusion of solute mass into the ice-rich permafrost (Figures 9a and 9b). Horizontally, the plume migrated approximately 225 m down gradient of the lagoon and 190 m upgradient of it. Compared to scenarios S1 and S2 with stable permafrost and shallow episodic groundwater flow (Figures 6 and 7), the conservative solute plume extent was greatly enhanced, and a much larger area around the lagoon was impacted. The solute mass from the lagoon however still did not reach the river in this simulation and was constrained at depth by the deepening permafrost table. It should be noted that nitrogen and phosphorus plumes were still much less persistent and are similar in extent to S1 and S2 despite the loss of permafrost.
In S4a, with an open talik underneath the lagoon, the conservative plume now extended through the talik to depths of 60 m below the surface of the domain (Figure 9f), but with the same horizontal extent as S3a.
Again, similar to S3a, the degrading and sorbing plumes were less persistent despite permafrost loss (Figures 9g and 9h). These results, along with those of scenarios S1 and S2 suggest that permafrost, if present and stable, can be an effective barrier to contaminant migration from the lagoon. As expected, thawing permafrost allows for the enhanced migration of solutes (Connolly et al., 2020), and the presence of taliks enhances the depth of solute migration (Johansson et al., 2015). However, the model shows that solute transformation processes like sorption and degradation, can significantly reduce solute migration despite permafrost degradation and increased groundwater flow and subsurface hydrologic connectivity.
Solute migration results for scenarios S3b and S4b with a high permeability in the upper 5 m of the domain (Figure 2c; Table 3) are shown in Figure 10. The increased permeability of this layer resulted in much more groundwater flow and thus more solute mass present in the system for conservative solutes (Figures 10b  and 10f) compared to the base case scenarios, but did not increased the plume's extent down gradient (Figures 9b and 9f). This was due to local groundwater discharge, consistent with the presence of wetlands and marshy areas at the field site (Section 3.1). The degrading nitrogen plumes are more persistent in these scenarios and extended approximately 150 m from the lagoon. In these scenarios, however, the conservative plumes still did not impact the river in the 100-year simulation period. It should be noted that the presence of a high permeability surficial layer in scenarios S3b and S4b (Figure 10) allowed the transport regime to become uncoupled from thaw dynamics, as the high conductivity kept most water flow within the near surface layer (upper 5 m), above the thawing permafrost table (Figures 10a and 10e). Scenarios S3 and S4 MOHAMMED ET AL.

10.1029/2020WR028771
13 of 20 highlight conditions under which permafrost presence is not the main driver of solute transport dynamics, and show that under certain conditions hydraulic (Figures 10b and 10f) and/or reaction properties (reactive plumes in Figures 9 and 10) may drive the fate and persistence of solutes rather than the presence of permafrost.
The S5 simulation scenario utilized a highly idealized model domain to purposefully impose hydrogeologic conditions that would maximize the rate of contaminant transport to the river. The model domain was reduced to 3 m in thickness to represent only the active layer, and no freeze-thaw was simulated so that the subsurface never froze and groundwater flow occurred year-long from the lagoon to the river. Under these conditions, breakthrough of the conservative solute at the river occurred at approximately 11 years (Figure 11a). For the degradable contaminant nitrogen, breakthrough never occurred (Figure 11b), and the plume reached a dynamic steady-state condition upgradient of the river. Figure 10c shows that the phosphorus (sorbing solute) plume migrated more than 600 m, but still had not reached the river after 500 years.

Contamination Risk to River
Based on simulations conducted in these scenarios (Table 3), the likelihood of contaminants associated with the wastewater lagoon at the study site reaching the river is low. Simulations suggest groundwater and surface water in the immediate vicinity (∼250 m) of the lagoon may be impacted by more conservative solutes (e.g., chloride). As the water table at the site is relatively close to surface where roots, fractures, and macropore networks exist, shallow preferential flow near the lagoon (which were not simulated) or thawing permafrost (which was simulated) may allow degradable solutes to affect this immediate vicinity as well. As the site is in a region of discontinuous permafrost with known permafrost degradation (Smith et al., 2009), results from scenarios S3 and S4 are more likely to represent reality at the field site. Evidence to support this include elevated concentration of solutes in the monitoring wells within 250 m of the lagoon (Supporting information S1). Areas of local groundwater discharge such as marshes and wetlands in the area may also be affected by persistent solutes (Figures 9 and 10).
It should be noted that the approach utilized here, in which only the saturated zone was modeled, does not consider the formation of a shallow perched aquifer and shallow subsurface flow which is the major runoff mechanism in many boreal forests in permafrost environments (Evans & Ge, 2017;Jafarov et al., 2018;Jan et al., 2020;Neilson et al., 2018). However, this pertains to snowmelt-driven runoff (surface and subsurface) and event flow generation, and not groundwater recharge to the permanently saturated zone and baseflow generation. Because the lagoon is located at a local topographic high, snowmelt from the surrounding landscape is unlikely to enter the lagoon, and thus neglecting this mechanism should not invalidate the main site-specific conclusions. However, we acknowledge the current model formulation and site conceptualization preclude consideration of these flow mechanisms, and future modeling should consider these processes to assess our assumptions.
Despite these limitations, all simulations, including "worst-case" scenarios, show that degradable solutes are likely naturally attenuated in this specific setting before discharge to the river. For example, in these simulations, nitrogen (degrading solute) or phosphorous (sorbing solute) plumes are unlikely to ever reach the river. This preliminary modeling suggests that if permafrost is present, this type of wastewater management approach could be a viable option for other similar small communities in the North. However, as scenarios here have shown, the performance of this type of system depends on a multitude of site-specific factors including topography, hydrogeology, soil properties, distance to surface water systems, climate, and permafrost (hydraulic and thermal) conditions. Also, systems that may perform well today may fail in the future due to climate change and concomitant permafrost thaw and activated groundwater flow systems.
MOHAMMED ET AL.

Model Limitations and Future Work
Key limitations of the model presented are: (a) it does not currently simulate vadose (unsaturated) zone processes; (b) it does not consider fluid density changes; and (c) it handles concentration-dependent freezing in a simplified manner that requires calibration. In addition to significantly reducing computational effort, these simplifications were considered acceptable in our example as the lagoon effluent concentrations were not expected to alter fluid density enough to affect results (e.g., the Raleigh stability criterion (Schincariol et al., 1997) was not exceeded). The addition of unsaturated zone flow and transport would allow for perched aquifer dynamics described above (Section 4.4) that were not represented here, and allow for more complex solute degradation and reaction processes to be accommodated in the model, such as thaw-induced mobilization of heavy metals and saturation-dependent organic matter decomposition from Arctic soils. Inclusion of variable-density flow and more physically based constitutive relations describing the effects of solutes on porewater freeze-thaw would enable the model to be applied to more diverse cold region environments such as coastal permafrost settings. Further research is underway to incorporate these processes and examine how they affect transport and reaction mechanisms in cold region hydrogeological systems. The difficulty of understanding the complex mechanisms governing contaminant/nutrient fate and transport in thawing permafrost environments necessitate more field and experimental data to validate emerging models like the one presented here.
There have already been significant advancements, albeit limited in number, in other methods of simulating solute transport and permafrost processes in groundwater models, such as particle tracking and Lagrangian travel-time based methods (Bosson et al., 2013;Frampton & Destouni, 2015). It should be noted that the simulations here neglected subsurface heterogeneities in flow, transport, and reaction properties, which travel-time methods have been able to efficiently capture, even with multi-component solute reaction and water-rock interaction processes (Cvetkovic et al., 2012;Malmström et al., 2004;Trinchero et al., 2020). The approach presented here which utilizes the advection-dispersion equation provides another useful tool for simulating reactive solute transport in permafrost groundwater systems. Furthermore, the methodology employed here where constitutive relationships were used to couple heat and solute transport equations and relate solute reactivity to temperature conditions and heat transfer, cannot currently be simulated with a particle-tracking approach. This approach can also be incorporated into other similar formulations of groundwater flow and advective-dispersive mass/energy transport models (e.g., McKenzie et al., 2007;Diersch, 2013;Dagenais et al., 2020). It should be noted that the coupling of multiple transport equations in this approach, is computationally burdensome compared to particle tracking methods. However, massively parallelized multi-physics hydrological models like ATS (Painter et al., 2016), PFLOTRAN-ice (Karra et al., 2014), and OpenFoam (Orgogozo et al., 2019), have robust numerical schemes and already have multiple heat transfer, solute transport, and multiphase flow routines incorporated, and thus also have the numerical structure to accommodate the coupled heat-solute transport processes presented here.

Summary and Conclusions
This study presented the development and application of a physically based numerical model for simulating solute transport in permafrost-affected groundwater systems. Simulations identified important coupled water-heat-solute transport mechanisms controlling the fate of groundwater contaminants in permafrost Figure 11. Scenario S5 simulation with no permafrost showing (a) conservative plume (10 mg/L isoline) reaching the river at approximately 11 years and breakthrough curve (inset); (b) total nitrogen (10 mg/L isoline) and (c) phosphorus (1 mg/L isoline) plumes at 500 years and never reaching the river. Note the nitrogen simulation reaches a dynamic steady state, whereas phosphorus was still migrating toward the river at 500 years. regions, including how pore ice limits pore-connectivity and associated advective-dispersive transport. A modified concentration-dependent SFC showed that when elevated solute concentrations decrease the freezing temperature of porewater, thawed-induced transport can occur and may be an important process that facilitates the growth of contaminant plumes in warm permafrost. In addition, seasonal freezing causes the temporary immobilization of water and solutes, which increases solute residence times in the active layer. These mobilization-immobilization cycles combined with different degradation properties lead to species-specific transport times and mass attenuation dynamics, which strongly influence solute degradation processes. Shallow flow and seasonal freezing in the active layer restrict connectivity of flow and transport pathways to regional discharge zones. This creates longer transport pathways and increased residence times, solute reactions, and thus contaminant degradation.
As the climate warms and permafrost degrades, transport pathways and hydrologic connectivity increase due to enhanced deep groundwater flow through taliks and enhanced supra-permafrost groundwater flow, which expectedly increases the transport of conservative solutes. However, this also increases temperature-dependent solute reactions and sorption to unfrozen aquifer material, which makes the fate of reactive solutes dependent on other hydrogeologic factors in addition to permafrost presence/absence. Further model development should consider surface water and vadose zone processes, as well as more complex subsurface heterogeneities in flow and transport parameters. Incorporating more physically based descriptions of solute effects on porewater freezing should also be investigated. Finally, more field-based research is required to validate predictions from these new models, and to develop site assessment and science-based guidelines specifically for waste management facilities operating in subarctic and Arctic regions.

Data Availability Statement
The SHAW model is publicly available free of charge at https://www.ars.usda.gov/pacific-west-area/boise-id/northwest-watershed-research-center/docs/shaw-model/. FlexPDE (version 7.19) is a commercial code that can be purchased and downloaded at https://www.pdesolutions.com/. FlexPDE scripts used in this study are available on a public repository (https://data.mendeley.com/datasets/924kr89tfg/1). Northwest Territories), Justin Hazenberg and Jamie Goddard (Municipal and Community Affairs, Government of Northwest Territories) and Heather Scott (Mackenzie Valley Land and Water Board) for their assistance and feedback. We also thank Dr. Julia Guimond (Dalhousie University) for assistance with running simulations. This work was funded by the Government of the Northwest Territories, and the Natural Sciences and Engineering Research Council of Canada through a NSERC Discovery grant. This research was conducted while A. Mohammed held a fellowship through the Ocean Frontier Institute, an award from the Canada First Research Excellence Fund. An earlier version of the manuscript greatly benefited from reviews of Dr. Ethan Coon and 2 anonymous reviewers.