Saltwater Intrusion Intensifies Coastal Permafrost Thaw

Surface effects of sea‐level rise (SLR) in permafrost regions are obvious where increasingly iceless seas erode and inundate coastlines. SLR also drives saltwater intrusion, but subsurface impacts on permafrost‐bound coastlines are unseen and unclear due to limited field data and the absence of models that include salinity‐dependent groundwater flow with solute exclusion and freeze‐thaw dynamics. Here, we develop a numerical model with the aforementioned processes to investigate climate change impacts on coastal permafrost. We find that SLR drives lateral permafrost thaw due to depressed freezing temperatures from saltwater intrusion, whereas warming drives top‐down thaw. Under high SLR and low warming scenarios, thaw driven by SLR exceeds warming‐driven thaw when normalized to the influenced surface area. Results highlight an overlooked feedback mechanism between SLR and permafrost thaw with potential implications for coastal infrastructure, ocean‐aquifer interactions, and carbon mobilization.

than 100%, saltwater can intrude into unfrozen pore space via diffusion and advection and can depress the porewater freezing point (Harrison & Osterkamp, 1978), triggering thaw ( Figure 1). In turn, thaw-driven increases in permeability can further facilitate SWI. This mechanism of SWI-driven thaw drives conversion of ice-saturated permafrost to ice-free permafrost, where ice-free permafrost is permafrost with high liquid water saturation (Figure 1). At the surface, saltwater inundation can cause density-driven convection of salt water into the underlying aquifer system (Baker & Osterkamp, 1988) and surface warming induced permafrost thaw (Romanovskii et al., 2005), processes that further threaten coastal permafrost. A transition to ice-free permafrost also decreases the sediment's mechanical strength (Hivon & Sego, 1995), increasing susceptibility of coastlines to erosion (Arenson et al., 2007) and land subsidence (Nelson et al., 2001), both of which can further drive SWI.
Predictive understanding of the mechanistic links and relative timing of SLR and permafrost thaw has been hampered by the difficulty of representing interrelated physics of water flow coupled to energy and mass transport processes in coastal permafrost environments. In particular, cold-regions multi-physics models for coastal zones need to represent the latent heat effects of freeze-thaw processes (Lamontagne-Hallé et al., 2020), as well as solute concentration impacts on density-dependent flow and water freezing temperatures. Models that incorporate water and energy transport with latent heat have provided powerful insights into the interactions between permafrost thaw and groundwater flow (Bense et al., 2009(Bense et al., , 2012Ge et al., 2011;Lamontagne-Hallé et al., 2018;Painter et al., 2016). For coastal marine environments, Frederick and Buffett (2014) developed a three-phase groundwater flow model that simulates water, heat and mass transport to assess subsea groundwater flow impacts on permafrost and gas hydrate stability (Frederick & Buffett, 2015). Current cryohydrogeological models, however, do not consider the effects of solute exclusion during subsurface freeze-thaw, which we hypothesize to be crucial for coastal permafrost due to salinity-effects on permafrost characteristics that govern groundwater flow (Zhou et al., 2020). The combination of recent advances in our understanding of freezing processes in saline soils (e.g., Zhou et al., 2020) and the development of cold-regions groundwater models for terrestrial settings (Grenier et al., 2018) has provided the theoretical foundation from which we develop a numerical model that more fully integrates transient subsurface physical processes relevant for coastal, cold regions hydrogeological dynamics.
In this study, we investigate mechanistic links between changing ocean, atmosphere, and cryosphere conditions in rapidly evolving Arctic coastal environments. To do this, we develop and apply a process-based cryohydrogeological model that enables investigation of the impacts of SLR, and the warming of land and sea on Arctic coastal permafrost. The following sections will (a) describe the numerical modeling approach, (b) assess the impacts of SLR and warming on coastal permafrost, and (c) discuss implications for coastal Figure 1. Schematic diagram of a coastal Arctic cross-section showing (a) initial dynamics before climate change, (b) transitionary dynamics as sea levels rise and land surface temperatures increase, and (c) future dynamics with prolonged climate change. White dashed line shows initial coastline with ocean on the right and land on the left. Processes 1-3 in the diagram drive processes 4-5 but are temporally concurrent. The vertical exaggeration is sixfold. communities and ecosystems. Although we focus on the interpretation of simulations to provide new insight into complex processes in archetypical coastal permafrost environments, the model development is an important step in itself to enable future research in these settings.

Numerical Model
Two-dimensional hydrological models were developed within the generic finite element modeling environment provided by FlexPDE (PDE Solutions Inc., 2020, http://www.pdesolutions.com), which was used to simultaneously solve coupled partial differential equations describing subsurface, variable-density fluid flow and heat transfer with salinity-dependent freeze-thaw. FlexPDE has previously been used in terrestrial studies simulating permafrost thaw, groundwater flow, and solute transport (Bense et al., 2012;Grenier et al., 2018;Mohammed et al., 2021).
The model assumed a saturated aquifer and density-dependent, Darcian fluid flow: where h f is equivalent freshwater hydraulic head (m), ρ f (kg/m 3 ) is fluid density, ρ 0 is freshwater density (kg/ m 3 ), g (m/s 2 ) is acceleration due to gravity, μ (kg/m s) is fluid viscosity, k is the intrinsic permeability of the porous medium (m 2 ), k rw (−) is a relative permeability exponential function that reduces permeability as pore ice is formed (McKenzie et al., 2007) (Table S2), ϕ is sediment porosity (m 3 /m 3 ), S w is water saturation (−), and β is aquifer compressibility (m·s 2 /kg).
Energy transport was modeled using the conduction-advection equation considering the effects of latent heat associated with the melting/freezing of pore ice (Bense et al., 2009): where T (°C) is subsurface temperature,  E q (m/s) is the Darcy flux, C a (J/m 3 ·°C) and κ e (W/m·°C) are volumetric heat capacity and effective thermal conductivity of the sediment/fluid/ice mixture, respectively, and θ w (m 3 /m 3 ) is liquid water content. C w (J/m 3 ·°C) and L i (J/m 3 ) are volumetric water heat capacity and volumetric latent heat of fusion, respectively.
Solute transport was simulated using the advection-dispersion equation (Bear, 1972): where C is the solute concentration (g/m 3 ), D (m 2 /s) is the hydrodynamic dispersion tensor and v i (m/s) is the linear porewater velocity ( v q i w   / ). Changes in water density due to temperature and salinity are described using a nonlinear equation of state (Nguyen et al., 2020): To account for the lower freezing point of saline water due to an increase in the osmotic potential, we incorporated a recently developed salinity dependent soil freezing curve, which describes the freezing point and phase transition between liquid water and ice based on both temperature and salinity (Zhou et al., 2020): where S wres is residual water saturation (−), T 0 (°C) is the freezing point of bulk water, K F (°C·L/mol) is the coefficient of freezing point depression that lowers the freezing point of saline porewater (Zhou et al., 2020), and a and b are fitting parameters (−) based on values from Zhou et al. (2020). Solute exclusion occurred as liquid water froze and ice was formed, which increased the salt concentration in the liquid water as ice released the salt molecules. The effect of solute exclusion was incorporated by calculating the effective ion molarity, c e (mol/L), in the equilibrium solution of soil (Zhou et al., 2020): where c (mol/L) is the ion molarity. All model parameters and relations are listed in Table S1.

Model Domain, Simulations, and Assessment
The model domain represented an archetypical, two-dimensional, vertically oriented coastal Arctic transect 1,250 m in length with 1,000 m on land (0-1,000 m) and 250 m offshore (1,000-1,250 m) ( Figure 2). Simulation results confirmed that domain limits were sufficient to represent the maximum extent of SWI and lateral permafrost thaw, while retaining reasonable computational runtimes. The surface boundary of the domain represented the water table, which is assumed to coincide with topography (e.g., Ge et al., 2011;Mc-Kenzie et al., 2007), with a uniform 1% slope, a typical slope based on a global water table slope assessment (Luijendijk et al., 2020). The model domain extended to a uniform depth of 200 m below sea level.
The hydraulic head along the top of the landward domain was fixed equal to surface elevation (h = elevation) to specify the water table position (Figure 2a), representing a topography-limited system where water table elevation does not respond to changing sea levels. Topography-limited conditions are found along coastlines Arctic-wide (Michael et al., 2013). The hydraulic head along the top of the seaward domain (i.e., inundated with seawater) was equal to sea level (h = sea level). The vertical sides and base were no-flow boundaries assuming a watershed divide along the landward boundary and ignoring potential effects of groundwater flow across the seaward and bottom boundaries (Figure 2a). Similarly, the vertical sides and base of the domain were closed to mass transport. The salt concentration boundary conditions along the top of the domain were spatially variable. On the landward domain, a specified water flux boundary at the surface was directionally dependent for concentration, with freshwater recharging the system and model-computed porewater concentrations discharging from the system. The surface of the seaward domain had a specified concentration boundary equal to that for seawater (35 g/L).
The domain surface was also assigned a specified temperature boundary. From 0 to 500 m, the temperature decreased linearly from 0.5 to −4°C, simulating the presence of a talik on the landward boundary of the domain (Figure 2a). The talik represented conditions under a lake, stream, river basin, or areas with a regional connection to recharge zones and enabled a regional sub-permafrost groundwater flow system to exist (Rowland et al., 2011). This deep groundwater flow has been suggested in other studies (Kane et al., 2013;Scheidegger & Bense, 2014) and is necessary to provide a counterforce to inhibit the saline front from sustained intrusion and associated domain-wide permafrost degradation, which would be inconsistent with present-day coastal permafrost presence. From 500 to 1,000 m, an initial temperature of −4°C was applied laterally to the domain top. On the seaward domain, the domain surface temperature increased linearly from −4 to 0°C from the coast (1,000 m) to the end of the domain (1,250 m) reflective of warmer submarine ground temperatures (Harris et al., 2017;Pedrazas et al., 2020;Romanovskii et al., 2005). The domain bottom had a constant geothermal heat flux of 0.05 W/m 2 , which is representative of median conditions across the Arctic (Davies, 2013), and the vertical sides were closed to heat exchange. Initial steady state was established by holding the top temperature boundary condition and sea level constant (h = 0) until the water flow, temperature, and ice-saturation stabilized. While other approaches have been used to 'spin-up' permafrost models for site-specific analyses (Frederick & Buffett, 2014), given the synthetic models used in this study, model simulation began from steady state. Steady-state conditions were established for three intrinsic permeability scenarios (k = 1 × 10 −13 m 2 , 1 × 10 −14 m 2 , 1 × 10 −15 m 2 ), representing unconsolidated and siliciclastic sedimentary deposits (Gleeson et al., 2011).
Prior to running fully transient simulations we tested the impact of solute exclusion on initial permafrost formation and distribution by running simulations with and without solute exclusion (i.e., with c e in Equation 5 or with c). The incorporation of solute exclusion strongly impacted the computed ice-saturated permafrost extent and saltwater distribution (Figures 2b and 2c). Without solute exclusion, ice-saturated permafrost formed under the saltwater boundary, and the saltwater wedge was flushed out. In contrast, with solute exclusion a permafrost wedge formed at the coast, creating initial ice-saturated permafrost distributions similar to those derived from recent geophysical measurements of such a system (Kasprzak, 2020). Additionally, the saltwater wedge was present and hypersaline conditions (i.e., salinity greater than seawater) were established as observed in high-latitude marine settings (Harris et al., 2017;Spirina et al., 2017). Given the evident importance of solute exclusion, all simulations discussed later include this process.
All transient simulations were run from steady-state initial conditions for a total of 120 years based on the time from the onset of rapid Arctic warming (i.e., 1980;Osterkamp, 2003Osterkamp, , 2005 to 2100 using an adaptive timestepping approach. Five sea-level change scenarios were considered including one sea-level fall (−0.005 m/year), one constant sea level, and three SLR (0.002, 0.004 and 0.008 m/year) scenarios based on 2100 SLR projections (Table S2) (James et al., 2014;Oppenheimer et al., 2019;Vousdoukas et al., 2017). Additionally, five surface warming scenarios were considered (Table S2): one without warming, one with inundation-induced warming from overlying seawater, and three rates of inundation plus land and ocean warming based on Representative Concentration Pathways (RCPs) emission scenarios 2.6, 4.5, and 8.5 (Collins et al., 2013). Sea level and temperature were changed at a constant temporal rate over the simulation. Seasonal variability in land surface temperatures was not considered due to the present study's focus on mechanisms of thaw and decadal changes in permafrost extent, which are assumed to exceed interannual variability impacts on permafrost (Bense et al., 2009). This assumption may influence the exact magnitude and timeline of change in either direction. However, by removing spatial and temporal complexity in land surface temperature, new process-based understanding is generated, unobscured by incorporation of shorter spatial-and temporal-scale processes. Sea-level changes were expressed by an increase or decrease in the seaward and near-shore domain specified head boundary condition. As sea levels fell, the specified head boundary transitioned from mean sea level to land elevation. With sea-level rise, the seaward boundary head increased and transgressed landward, converting land elevation boundaries to sea-level boundaries. Similar changes were implemented for the concentration boundary condition, with the specified flux moving seaward during sea-level fall, and the specified salt concentration (35 g/L) moving landward during sea-level rise.
The projected global average air temperature increases for RCP 2.6, 4.5, and 8.5 scenarios were used to represent Arctic land temperature increases. This is conservative because recent Arctic air temperature increases are nearly double the global average (Collins et al., 2013). We recognize that this assumption is a limiting simplification of the model, as air temperature is only a proxy for land surface temperature and offsets between air and land surface temperatures are spatially and temporally variable (Mann & Schmidt, 2003). Numerous factors influence land surface temperatures including snow cover, surface water reservoirs, and vegetation (Jafarov et al., 2018;Rushlow et al., 2020;Young et al., 2020). Given our focus on decadal trends with warming and sea-level change using archetypal models, we chose readily available global average air temperature trends (Collins et al., 2013) to enable mechanistic assessment of the processes of interest. Changes in sea level were accompanied by changes in specified temperature due to inundation-induced warming. Land and sea warming were included in combination with inundation-induced changes in land temperature (Table S2). All sea-level change scenarios were run for all warming scenarios for a total of 25 simulations for each k scenario. Results are reported for k = 1 × 10 −13 m 2 unless indicated otherwise.
The model output was used to assess ice-saturated and total (<0°C) permafrost volume per meter coastline (m 3 /m) at every timestep. Results were reported for a subset of the model domain from 750 to 1,250 m (Figure 2). Permafrost was considered ice-saturated if the ice-saturation exceeded 30% of pore space. An ice-saturation threshold of 30% is conservative with respect to hydrogeological activation (Wellman et al., 2013). Given that warming drives permafrost thaw domain-wide and SLR impacts are confined to the coastal zone, we calculated the volume of permafrost change relative to the area of influence (m 3 /m 2 ) to compare thaw volumes between drivers. The relative volume of ice-saturated permafrost loss per area of inundation due to SLR (herein referred to as area-normalized thaw) was calculated by dividing the difference in permafrost loss per m coastline between warming with SLR (m 3 /m) and warming only (m 3 /m) by the surface area inundated by SLR per m coastline (m 2 /m). For example, the difference was divided by 24 m 2 /m (= 0.002 m 2 /m·years × 120 years), 48 m 2 /m, and 96 m 2 /m for low, moderate, and high rates of SLR, respectively. For warming only scenarios, loss per m 2 influenced was calculated by dividing the ice-saturated permafrost loss by 250 m 2 /m (i.e., 750-1,000 m along domain length × 1 m domain width). These simulations did not consider erosion impacts on permafrost extent (Fritz et al., 2017).

Sea-Level Rise and Warming Affects Coastal Permafrost
SLR reduced the coastal ice-saturated permafrost extent due to the landward migration of the subsurface freshwater-saltwater interface and resulting depression in freezing temperature (Figures 3 and 4a). SLR alone decreased ice-saturated permafrost extent by as much as 1,209 m 3 per m of coastline or an area-normalized loss of 12.6 m 3 /m 2 (Figures 3 and 4a).
The combination of SLR and land-and sea-warming decreased ice-saturated permafrost extent for all warming and sea-level change scenarios relative to pre-warming (i.e., 1980) conditions. In general, the magnitude of loss increased when rates of SLR and warming were higher (Figure 3). Under warming from the RCP 8.5 emission scenario and 0.008 m/yr of SLR, the projected ice-saturated permafrost loss by 2100 (i.e., 120 years) was 2616 m 3 per m coastline which translates to a 9.5% loss (Figure 3b). In contrast to the lateral thaw driven by SWI, warming resulted in top-down thaw (Figures 3, S1, and S2). The combination of lateral and vertical permafrost loss contributed to extensive coastal ice-saturated permafrost loss in the upper portion of the aquifer. Significant permafrost loss was not observed in the lower part of the aquifer.
Increasing surface temperature generally had a greater impact on permafrost loss than SLR because warming was domain-wide rather than confined to the coastal zone. However, when considering the area influenced by both processes at play (i.e., saltwater inundation and warming), the area-normalized permafrost loss due to SLR alone exceeded loss due to warming for select scenarios. For example, high rates of SLR (i.e., 0.008 m/yr) resulted in a greater area-normalized loss than warming for low and moderate warming scenarios ( Figure 4b). In contrast, for RCP 4.5 scenarios with low and moderate rates of SLR and all RCP 8.5 scenarios, the area-normalized permafrost loss by warming exceeded loss by SLR (Figure 4b). In RCP 2.6 scenarios, loss due to warming exceeded loss due to low rates of SLR, but not moderate and high SLR rates. These results suggest that SLR-driven subsurface SWI, an unseen and previously unexplored driver of coastal permafrost thaw, can exert a critical control on coastal permafrost degradation and can result in more loss of ice-saturated permafrost along the coast than terrestrial surface warming.
Total permafrost volume also decreased with SLR and warming combined (Figures 4c and S2, Table S3). However, warming had a greater impact on total permafrost than SLR alone. When sea-level change was ignored, by 2100 surface warming alone decreased permafrost volume by 6.2%, 12.2%, and 25.3% for RCP 2.6, 4.5, and 8.5 scenarios, respectively. The combination of SLR and warming decreased total permafrost by up to 30.5% for the highest SLR (0.008 m/yr) and warming (RCP 8.5) scenario as a consequence of top-down thaw and landward ocean transgression.

Intrinsic Permeability Affects Ice-Saturated Permafrost Loss
Ice-saturated permafrost loss was reduced in the models for lower unfrozen aquifer permeability (k) (Figure 4d) but exhibited similar patterns with sea-level change and warming for most scenarios. For example, for the same SLR (0.004 m/yr) and warming (RCP 2.6) scenarios, the ice-saturated permafrost loss was 380, 249, and 190 m 3 per m coastline for k values of 1 × 10 −13 m 2 , 1 × 10 −14 m 2 , 1 × 10 −15 m 2 , respectively. However, under the highest warming and lower k scenarios, ice-saturated permafrost loss decreased with SLR due to reduced SWI and thermal insulation of overlying seawater. Colors indicate sea-level rise (SLR) and warming scenario, where red is permafrost extent without warming or SLR; green is permafrost extent with warming only; and blue is permafrost extent with warming and 0.008 m/yr of SLR. Filled colors overlap (e.g., all of the blue zone is within the red zone). Initial coastline is at distance of 1,000 m.

Discussion and Conclusions
Amplified high-latitude climate change has been observed for decades and is projected to continue (Collins et al., 2013). Understanding of the effects of climate change on Arctic coastal hydrogeological environments lags that in temperate and tropical climate zones, largely due to limitations in field access and modeling capabilities. The model developed and applied in this study is a significant and necessary step to better understand the current rapid evolution of coastal Arctic groundwater systems. Model results strongly suggest that SLR-driven SWI and consequent changes in subsurface freezing temperatures reduce coastal ice-saturated permafrost extent. SLR has the potential to laterally thaw ice-saturated permafrost and, when combined with warming-driven thaw, can result in extensive coastal permafrost loss.
Ice-saturated permafrost loss has implications for coastal ecosystems and communities. For example, permafrost thaw decreases the strength of the sediment (Bull et al., 2020;Jorgenson et al., 2006). Thus extensive lateral thaw can enhance erosion rates and decrease coastal land stability with implications for the resilience of coastal communities and infrastructure they rely upon (Hjort et al., 2018). Furthermore, permafrost thaw activates groundwater flow systems and may influence aquifer-ocean exchange and water quality.
Permafrost thaw also has consequences for carbon mobilization and global climate. Upon thaw, large quantities of biolabile organic carbon may be released to the atmosphere (Schuur et al., 2015;Vonk et al., 2013;Walvoord et al., 2019), resulting in a positive feedback loop between permafrost thaw, atmospheric carbon, and global warming (Schuur et al., 2015). This study shows that incorporation of advection, salinity-dependent freeze-thaw, and solute exclusion will enable integrated assessments of coastal Arctic permafrost and hydrogeological processes and enhance understanding of future coastal conditions. Future efforts will explore implications of intensified thaw on carbon transport and groundwater discharge.