Sea-level rise and warming mediate coastal groundwater discharge in the Arctic

Groundwater discharge is an important mechanism through which fresh water and associated solutes are delivered to the ocean. Permafrost environments have traditionally been considered hydrogeologically inactive, yet with accelerated climate change and permafrost thaw, groundwater flow paths are activating and opening subsurface connections to the coastal zone. While warming has the potential to increase land-sea connectivity, sea-level change has the potential to alter land-sea hydraulic gradients and enhance coastal permafrost thaw, resulting in a complex interplay that will govern future groundwater discharge dynamics along Arctic coastlines. Here, we use a recently developed permafrost hydrological model that simulates variable-density groundwater flow and salinity-dependent freeze-thaw to investigate the impacts of sea-level change and land and ocean warming on the magnitude, spatial distribution, and salinity of coastal groundwater discharge. Results project both an increase and decrease in discharge with climate change depending on the rate of warming and sea-level change. Under high warming and low sea-level rise scenarios, results show up to a 58% increase in coastal groundwater discharge by 2100 due to the formation of a supra-permafrost aquifer that enhances freshwater delivery to the coastal zone. With higher rates of sea-level rise, the increase in discharge due to warming is reduced to 21% as sea-level rise decreased land-sea hydraulic gradients. Under lower warming scenarios for which supra-permafrost groundwater flow was not established, discharge decreased by up to 26% between 1980 and 2100 for high sea-level rise scenarios and increased only 8% under low sea-level rise scenarios. Thus, regions with higher warming rates and lower rates of sea-level change (e.g. northern Nunavut, Canada) will experience a greater increase in discharge than regions with lower warming rates and higher rates of sea-level change. The magnitude, location and salinity of discharge have important implications for ecosystem function, water quality, and carbon dynamics in coastal zones.


Introduction
Groundwater discharge is an important mechanism through which carbon-, nutrient-, trace-metal-, and contaminant-rich water is delivered to the coastal ocean (Luijendijk et al 2020, Mayfield et al 2020. This groundwater discharge to the coastal ocean occurs due to both oceanic and terrestrial drivers including land-sea hydraulic gradients, tidal pumping, and density-driven convection (Taniguchi et al 2019). A portion of discharged groundwater, known as circulated submarine groundwater discharge (SGD), originates from and discharges back to the ocean. Groundwater discharge can also originate on land and discharge to the coastal zone above or below mean sea level. This land-derived discharge is the focus of this study and is herein referred to as coastal groundwater discharge (CGD) (Luijendijk et al 2020).
While often fresh, CGD can vary in salinity depending on the groundwater flow path (e.g. transport through the subsurface freshwater-saltwater mixing zone) (figure 1). In recent decades, CGD has been a topic of great interest due to its implications for the health and functioning of coastal ecosystems (Moore 2010, Taniguchi et al 2019, Guimond and Tamborski 2021, Santos et al 2021. However, despite comprising 34% of the global coastline (Lantuit et al 2012), coastal permafrost regions have been largely overlooked in quantitative and qualitative CGD assessments. Recent work on the hydrologic transformations of inland permafrost landscapes (Walvoord and Kurylyk 2016) raises critical but unanswered questions regarding the future of high-latitude coastal hydrology in a changing climate.
Climate change has the potential to alter the magnitude, salinity, and location of groundwater discharge in coastal settings due to changes in sea level, water table elevation (i.e. through recharge and groundwater table response to sea-level rise (SLR)), and tidal range (Michael et al 2013, Haigh et al 2020. For example, Michael et al (2013) projected decreased fresh and saline groundwater discharge with SLR in topography-limited systems due to decreased landsea hydraulic gradients, and, depending on the recharge rate, landward migration of discharge zones with SLR. Discharge responses to climate change along Arctic coastlines are even more complex due to the combined effects of SLR and warming-driven permafrost thaw, which transforms the subsurface permeability distribution. Arctic atmospheric warming is double the global average, resulting in rapid and widespread permafrost thaw (Collins et al 2013), with an additional 50% loss in global permafrost area projected by 2080 under the representative concentration pathway (RCP) 4.5 emissions scenario (Overland et al 2019). Permafrost thaw activates or enhances groundwater flow through increases in active layer thickness and formation of perennially thawed zones (taliks), resulting in increased groundwater discharge to rivers (Smith et al 2007, Walvoord and Striegl 2007, Bense et al 2009, Jacques and Sauchyn 2009, Ge et al 2011, Duan et al 2017, Lamontagne-Hallé et al 2018, Walvoord et al 2019. Based on these findings, we hypothesize that groundwater discharge dynamics are also changing along northern coastlines; however, these dynamics are additionally influenced by processes (e.g. SLR and density-dependent flow) that do not govern in inland permafrost landscapes.
In the coastal zone, SLR has the potential to produce offsetting effects on coastal discharge. SLR increases the depth of the water column (i.e. pressure head) on the ocean boundary, which, in topographylimited systems, decreases the hydraulic gradient and, in turn, reduces groundwater fluxes between terrestrial water reservoirs and the ocean (Michael et al 2013). Additionally, a recent study using the model discussed in this work identified links between SLR and ice-saturated permafrost extent, as saltwater intrusion into unfrozen pore space decreases the permafrost freezing temperature, triggering thaw (Guimond et al 2021a). Thus, along permafrostbound coastlines SLR has the potential to both decrease hydraulic gradients (decreasing CGD) and enhance lateral permafrost thaw, thereby activating groundwater flow (increasing CGD). Simultaneously, top-down permafrost thaw enhances aquifer-ocean connectivity near the ground surface (increasing CGD). These competing processes make evaluation and projections of CGD rates and drivers challenging and so far, unexplored.
In this study, we use a recently developed permafrost groundwater model to assess climate change impacts on the magnitude, location, and salinity of CGD. Specifically, we conduct numerical experiments of a generalized land-ocean transect to address how variable rates of SLR/fall and surface warming impact the magnitude and location of CGD under multiple permafrost temperatures and aquifer permeability scenarios. Findings elucidate trends in the magnitude and patterns of groundwater delivered to the coast over the next century.

Numerical model
Two-dimensional groundwater flow models of a generalized coastal Arctic transect were developed using FlexPDE (PDE Solutions Inc., 2020, www. pdesolutions.com), a finite element software that solves the coupled partial differential equations describing variable-density groundwater flow and heat transfer with salinity-dependent freeze-thaw. The model equations and boundary conditions are detailed in Guimond et al (2021a), which only focused on permafrost thaw and saltwater intrusion, and will briefly be discussed here. The present study incorporates new scenarios and output analyses of groundwater discharge dynamics in space and time.
The model simulated a water-saturated domain with density-dependent, Darcian fluid flow and permeability that decreases exponentially with ice content (Bense et al 2009). Energy transport was modeled using the conduction-advection equation with latent heat effects associated with the melting/freezing of pore ice (Bense et al 2009). The advectiondispersion equation, modified for the presence of pore-ice, was used to simulate solute transport (Mohammed et al 2021). Changes in water density due to temperature and salinity were described using a nonlinear equation of state (Nguyen et al 2020). A newly developed salinity-dependent soil freezing curve (SFC) (Zhou et al 2020) was also included to model the phase transition between liquid water and ice considering the impact of salinity on the freezing point of porewater. This salinitydependent SFC also incorporates solute exclusion, Figure 1. Diagram showing (a) cross-sectional model domain with initial ice-saturated permafrost distribution, 0 • C isotherm, and freshwater-saltwater boundary (interface between blue shades), (b) initial groundwater discharge zones before SLR and warming, and (c) discharge zones following warming. A-E designations in diagram (a) are discussed in text. Diagrams (b), (c) delineate circulated SGD and CGD and highlight how CGD may be saline due to transport through the mixing zone. The grey, dashed line is the surface coastline projected vertically into the subsurface which moves in position between (b) and (c) due to SLR.
which occurs as liquid water freezes and expels salt molecules, increasing the salt concentration in the surrounding liquid water, and lowering the porewater freezing temperature. Table S1 (available online at stacks.iop.org/ERL/17/045027/mmedia) lists all model parameters and relations.

Model domain, simulations, and assessment
The model domain was 1250 m long with 1000 m onshore and 250 m offshore and extended to a depth of 200 m ( figure 1(a)). Simulation results confirmed that the offshore domain limits were sufficient to capture CGD without interference with the vertical boundary. The domain surface represented the water table with a uniform 1% slope across the domain surface, which maintained saturated conditions and is typical of much of the Arctic (Luijendijk et al 2020).
Boundary conditions were established for hydraulic head, mass transport, and temperature across each domain edge (i.e. surface, bottom, and sides). Specified head along the top terrestrial boundary (0-1000 m) was equivalent to the domain's surface elevation (h = elevation) (figure 1(a): (A)-(C)). The configuration enables simulation of a topography-limited system where the water table remains constant with changing sea level, a scenario proposed along many Arctic coastlines (Michael et al 2013). The surface mass transport along the terrestrial portion of the domain was directionally dependent: freshwater (0 g L −1 ) recharged the domain, while discharge concentrations were determined by modelcomputed concentrations. The seaward portion of the domain surface (initially 1000-1250 m) had a specified head and specified concentration boundary equal to sea level (h = sea level; initially h = 0) and seawater (35 g L −1 ), respectively (figure 1(a): (C)-(D)).
The temperature boundary conditions across the domain surface were spatially more complex. Between 0 m and 500 m (figure 1(a): (A)-(B)), the temperature decreased linearly from 0.5 • C to −4 • C, establishing a vertical talik along the inland side. The talik represented conditions under a lake, river or regional recharge area. The presence of a vertical talik was necessary to establish a regional groundwater flow system below the permafrost, as has been previously documented (Kane et al 2013, Scheidegger andBense 2014), and to prevent unrealistic domain-wide seawater encroachment and permafrost degradation. Sub-permafrost groundwater flow in the seaward direction counteracted saltwater intrusion, establishing present-day coastal permafrost presence. In reality the distance from ocean to vertical talik can range from hundreds of meters to hundreds of kilometers. However, given our focus on understanding coastal discharge responses under a range of climate change scenarios, we have avoided site-specific complexity and limited domain size due to the computational requirements (see section 4.3 for further discussion).
The initial surface temperature from 500 m to 1000 m ( Simulations were run in two steps. First, steady state was established by holding head, mass, and temperature conditions along the top of the domain constant until the water flow, temperature, and icesaturation fields stabilized (figure S1). Because these are not site-specific models with long-term climate records, we did not 'spin up' models with reconstructed climates as in selected other studies (e.g. Frederick and Buffett 2014). Steady-state conditions were established for three unfrozen intrinsic permeability scenarios (k = 1 × 10 −13 m 2 , 1 × 10 −14 m 2 , 1 × 10 −15 m 2 ) and two surface temperature scenarios (initial temperature = − 4 • C, −2 • C). Permeability scenarios represent fine-grained unconsolidated and siliciclastic sedimentary deposits (Gleeson et al 2011), characteristic of a majority of Arctic coastlines (Lantuit et al 2012, Overduin et al 2014. From steady state, simulations were run for 120 years  under five sea-level change and six warming scenarios (table 1). Sea-level change scenarios included rates of −5 mm yr −1 , 2 mm yr −1 , 4 mm yr −1 , and 8 mm yr −1 based on Arctic SLR projections for 2100 (James et al 2014, Vousdoukas et al 2017, Oppenheimer et al 2019 and a control scenario without sea-level change. Scenarios are referred to as no SLR, sea-level fall, and low, moderate, and high rates of SLR. Four warming scenarios include global average air temperature warming rates of 0.012 • C yr −1 , 0.024 • C yr −1 , and 0.048 • C yr −1 , respectively for the RCP emission scenarios 2.6, 4.5, and 8.5 (Collins et al 2013), as well as a control scenario without warming. The use of global average air temperature rate of change is conservative based on recent soil temperature projections for the Arctic region (Soong et al 2020). Thus, to capture the range of potential warming rates for permafrost soils, two additional warming scenarios were run with RCP 4.5 and 8.5 rates based on warming for shallow permafrost soils (0.033 • C yr −1 , and 0.062 • C yr −1 ; Soong et al 2020). The warming scenarios are herein referred to according to the associated RCP scenario, with the two RCP 4.5 and 8.5 scenarios delineated by low (based on global air temperature projections) and high (based on permafrost soil temperature projections) designations (table 1). Sea-level changes were accompanied by landward shifts in specified temperature due to seawater inundation. The slope of the offshore temperature boundary remained unchanged but shifted landward with SLR. Land and sea warming were included in combination with inundation-induced changes in land temperature. Together these changing hydraulic and thermal boundary conditions resulted in 30 runs for each intrinsic permeability and initial temperature scenario (total 120 runs).
Model output was assessed for fresh, saline, and total CGD (i.e. terrestrially derived discharge), as well as the location of total groundwater discharge (TGD, i.e. CGD plus circulated SGD, figures 1(b) and (c)) in the coastal zone. To calculate the discharge, the groundwater discharge flux and salinity at every node along the domain surface were output every two years for each simulation. The discharge was then integrated along segments of the domain to assess total discharge at a given time. TGD was assessed by integrating between 750 m and 1250 m along the domain length, capturing all CGD and circulated SGD. CGD was calculated by subtracting the circulated SGD (equivalent to offshore recharge) from the TGD. CGD was partitioned based on the nodal salinity along the domain surface, where fresh groundwater discharge was all discharge with less than 10 g L −1 of salt, and saline discharge was all discharge with a salt concentration of 10 g L −1 or greater (Post et al 2013). Results are reported for a subset of the model domain from 750 m to 1250 m, and, unless otherwise Table 1. Sea-level change, and land and ocean warming rates included in this study. All sea-level change scenarios were run with each warming scenario for a total of 30 simulations for each intrinsic permeability and initial surface temperature scenario. Warming rates are associated with the RCP 2.6, 4.5 (low), and 8.5 (low) scenarios (Collins et al 2013). RCP 4.5 (high) and 8.5 (high) are associated with projected permafrost soil warming scenarios (Soong et al 2020).

Climate change impacts on coastal groundwater discharge
Climate change impacted the quantity of fresh, saline and total CGD with the direction and magnitude of change dependent on the interplay between warming rates and SLR. Fresh CGD increased under higher warming rates (RCP 4.5 (high), RCP 8.5 (high and low); figure 2(a)). In 1980, fresh CGD accounted for less than 1% of CGD, but by 2100, fresh CGD accounted for up to 69% of discharge, a two-order-ofmagnitude increase. This increase is attributed to lateral talik formation due to top-down permafrost thaw under the highest warming scenarios (RCP 4.5 (high) and RCP 8.5 (high and low); figures 3(c) and (d)), opening a supra-permafrost conduit for groundwater flow and a subsurface connection to the coastal ocean. With a warming rate of 0.048 • C yr −1 (RCP 8.5 (low)), talik formation occurred 32 years and 74 years after 1980 (i.e. 2012 and 2054) for −2 • C and −4 • C initial surface temperature scenarios, respectively. The impacts of SLR on fresh CGD were more complex. For the RCP 2.6, 4.5 (low), and no warming scenarios, fresh discharge increased with SLR. For example, without warming, fresh CGD increased from 1.8 × 10 −5 m 3 d m −1 without SLR to 3.3 × 10 −3 m 3 d m −1 under the highest SLR scenario. Given that a lateral talik did not form in scenarios with low or no warming (figures 3(a) and (b)), this increase can be attributed, at least in part, to SLR enhanced lateral permafrost thaw (Guimond et al 2021a). Scenarios with higher warming rates (RCP 4.5 (high), RCP 8.5 (high and low)) showed the greatest fresh CGD with moderate SLR rates, but fresh discharge decreased slightly at high SLR rates (figure 2(a)), suggesting that the increased sea level for high SLR rates decreased the land-sea hydraulic gradient and overcame the effects that amplified CGD at lower rates of SLR.
Saline CGD was minimally impacted by warming alone, with the most notable impacts from the highest warming scenario which showed a decline in saline CGD ( figure 2(b)). This decrease is a result of enhanced groundwater discharge routed through the supra-permafrost aquifer in a more hydrologically active shallow flow system. Unsurprisingly, SLR decreased saline CGD due to the increased water head at the ocean boundary ( figure 2(b)).
In total, CGD discharge increased between 1980 and 2100 with higher warming and lower SLR rates. For example, CGD increased between 1980 and 2100 for all RCP 8.5 scenarios (low and high) but decreased for three of five RCP 2.6 scenarios (figure 2(c)). SLR had a contrasting effect on the change in discharge between 1980 and 2100, as observed by scenarios with higher SLR showing the largest decrease or smallest increase in CGD (figure 2(c)). Scenarios with sealevel fall resulted in increased discharge between 1980 and 2100 for all warming scenarios.
The interplay between warming and SLR mediated projected discharge rates with results described here for 2100. Warming from the RCP 2.6 and 4.5 (low) emission scenarios only minimally influenced CGD ( figure 2(a)). RCP 4.5 (high) and 8.5 (low and high) scenarios had the highest CGD due to the formation of a lateral talik that enabled supra-permafrost groundwater flow and connectivity between groundwater reservoirs and the coastal ocean ( figure 2(c)). The talik depth and consequent discharge was greatest for the RCP 8.5 (high) scenario. For a given warming rate, CGD decreased with SLR for all warming scenarios due to the flattening of the hydraulic gradient. The impact of SLR on CGD was greatest for the RCP 4.5 (high) and RCP 8.5 (low) warming scenarios where lateral talik formation occurred, as observed by a greater slope between the no SLR and highest (8 mm yr −1 ) SLR scenarios (figure 2(c); RCP 4.5 (high) slope = 9 × 10 −4 ; RCP 2.6 slope = 5 × 10 −4 ). However, with continued talik formation and overall greater discharge, the impact of SLR decreased as indicated by a decrease in the slope between no SLR and the highest SLR scenarios for RCP 8.5 (high) (slope = 6 × 10 −4 ) compared to RCP 4.5 (high) (slope = 9 × 10 −4 ) and RCP 8.5 (low) (slope = 8 × 10 −4 ). Figure 2 shows results for scenarios with an initial temperature of −4 • C. Scenarios with an initial temperature of −2 • C exhibited similar trends in CGD with SLR and warming, but the magnitude of discharge was greater ( Figure S2(c)). For example, the average CGD at 2100 for all runs with −2 • C initial surface temperature was 0.058 m 3 d m −1 , whereas the average for all runs with −4 • C initial surface temperature was 0.019 m 3 d m −1 . Unsurprisingly, a decrease in permeability had the opposite effect and resulted in order-of-magnitude or greater decrease in CGD (figures S2(a) and (b)). For example, the average discharge at 2100 was 1.6 × 10 −3 m 3 d m −1 and 1.2 × 10 −4 m 3 d m −1 for unfrozen k = 1 × 10 −14 m 2 and k = 1 × 10 −15 m 2 , respectively.

Changes in the location of CGD with climate change
Climate change also impacted the groundwater discharge location in the coastal zone through changes in hydrology and permafrost distribution (figure 4). While warming drives top-down thaw, SLR-driven saltwater intrusion drives the lateral retreat of icesaturated permafrost (i.e. permafrost with over 30% of pore space as ice) in the coastal zone (Guimond et al 2021a) (figure 4(c)). With warming, discharge shifted landward, even without SLR, due to suprapermafrost groundwater flow, resulting in discharge landward of mean sea level (figures 3 and 4(a)). Discharge shifted 41 m landward between the no warming and highest warming scenario. The lateral retreat of permafrost with SLR also expanded the groundwater flow field along the coast. With an increase in SLR, the location of maximum discharge transgressed Figure 3. Cross-section of ice-saturated permafrost extent at 1980 and 2100 and freshwater-saltwater interface at 2100 with groundwater flow paths with low (2 mm yr −1 ) and high (8 mm yr −1 ) SLR, RCP 2.6 and RCP 8.5 (high) warming, and k = 1 × 10 −13 m 2 scenarios. D = climate warming-driven groundwater discharge; ∆Q = change in discharge due to decreased land-sea gradients. landward due to ice-saturated permafrost thaw and hydrogeologic activation (figures 3 and 4(b)).

Changes in the quantity of CGD due to climate change
These simulations suggest that groundwater discharge will undergo pronounced changes in the coastal zone in topography-limited systems, with the magnitude and direction of change dependent on local warming and sea-level change. Results indicate that coastlines with the highest warming and lowest SLR projections will experience the greatest increase in CGD due to hydrogeological activation from topdown permafrost thaw and higher land-sea hydraulic gradients ( figure 3(c)). This includes regions such as Canadian Arctic coastlines which are experiencing rapid warming and sea-level fall due to isostatic uplift (Collins et al 2013, James et al 2014, Oppenheimer et al 2019. In contrast, these simulations also suggest that areas with lower warming and higher SLR rates, such as the Alaskan northern coastline, may experience a decrease in discharge due to shallower zones of supra-permafrost groundwater flow and decreased land-sea hydraulic gradients (figure 3(b)) (Collins et al 2013, Oppenheimer et al 2019. However, much more research is needed to determine how local factors such as sediment type, heterogeneity, and erosion impact CGD projections. This study is an important first step towards better constraining coastal high-latitude carbon dynamics, as changing discharge patterns will impact the lateral delivery of thaw-mobilized carbon to the coastal ocean. Sourced from the supra-permafrost aquifer, fresh CGD may contain significant concentrations of newly mobilized, terrestrially-derived carbon and nitrogen (Connolly et al 2020), mercury (Schuster et al 2018), or other sequestered constituents (Legendre et al 2014). However, the fate of this transformed and translocated carbon and nutrients once delivered to the coastal zone is uncertain. Some carbon may be released to the atmosphere as greenhouse gas, contributing to the positive permafrost-carbon feedback with global warming (Schuur et al 2015); some of the carbon may also  Guimond et al (2021a). Dashed lines in (b) indicate mean sea level location at 2100 for the corresponding sea-level change scenario. The original shoreline location at 1000 m is indicated by the black, dashed line. We hypothesize that the coastal recharge in (b) for the 8 mm yr −1 scenario is due to steep pore ice and salinity gradients resulting in a coastal circulation zone, but further research on coastal groundwater flow paths is needed. be consumed in the coastal ocean, fueling phytoplankton biomass growth and potential eutrophication (Rodellas et al 2014, Adolf et al 2019; and some may be buried in sediments (Jong et al 2020).

Permafrost impacts on groundwater discharge location
The landward shift in discharge location with climate change could trigger changes in Arctic ecosystem ecology and biogeochemistry. CGD delivers solutes, including carbon, nutrients, and metals, directly to the ocean, where they enter the food web and impact coastal water quality (Taniguchi et al 2019, Santos et al 2021. CGD that is delivered onshore meets the atmosphere and/or surface water reservoirs (i.e. tidal rivers, channels, ponds, wetlands) and is potentially transpired or transformed before/during mobilization to the coastal ocean. Even if delivered to the ocean via overland flow pathways, the longer flow paths and potentially aerobic surface water environments change the biogeochemistry of the ocean-reaching water.
Similarly, lateral talik formation and CGD landward of mean sea level may support the development or maintenance of coastal freshwater wetlands. For example, Lamontagne-Hallé et al (2018) found that shallow lateral taliks in inland settings support wetland-like conditions along rivers, but with increasing thaw, discharge shifts towards the river. Thus, landward discharge and wetland-like conditions are temporary and indicative of a system in transition. Our models which include SLR and density-dependent flow dynamics support similar findings but for the coastal zone, with increased freshwater discharge to the coastal zone with time ( figure 4).

Uncertainties and limitations
These results provide the first insight into climate change-driven changes in high-latitude coastal groundwater discharge. The generalized approach enables process-based understanding across a range of warming and SLR scenarios but requires simplifications due to the dearth of field data and computational demand of simulating variable-density groundwater flow with salinity-dependent freezethaw. First, the simulations assumed a saturated domain and constant terrestrial groundwater table (i.e. topography limited) that is unresponsive to SLR induced changes. We assume an increase in water table elevation would intercept the land surface and discharge via overland flow pathways. This approach is computationally more stable and less demanding than simulating unsaturated groundwater flow. However, in some areas of the Arctic, such as the rechargelimited regions of the Arctic Archipelago (Michael et al 2013), groundwater table rise in response to SLR will not intercept the land surface, maintaining landsea hydraulic gradients and reducing saltwater intrusion (Befus et al 2020). Given that our model domain falls within 1 km of the ocean where there is minimal difference between topography-and recharge-limited systems due to shallow water tables (Luijendijk et al 2020), we do not expect our topography-limited assumption to significantly impact our results. Conversely, coastal watersheds are often laterally more extensive than 1 km, which would increase discharge if recharge scaled proportionally and impact the tipping point between climate warming-driven groundwater discharge and the change in discharge due to decreased land-sea gradients. Future studies are needed to understand discharge changes with different geometric variations.
The water table boundary also mediates the magnitude and location of recharge to the model domain. In general, modeled recharge ranges between 0 and 100 mm yr −1 , in line with Arctic estimates (Döll and Fiedler 2008). Rather than distributed across the domain, recharge was predominantly confined to the landward edge of the domain, an artifact of the water table boundary and vertical talik. Future studies and model developments that integrate a recharge and seepage boundary will be critical to refining the magnitude of change, as well as incorporating dynamic unsaturated zone processes and seasonally varying land surface temperatures guided by observational CGD data across variable permafrost, climate, and hydrogeological regimes.
Lastly, in this study we represent an isotropic, homogeneous model of a generalized land-sea transect to assess changes with SLR and warming. Coastal aquifers can have heterogeneity that transcends spatial scales from confining units, sediment type, and macropores, which impact permeability, saltwater intrusion, groundwater flow paths, and ultimately CGD. This study incorporates conservative permeability scenarios (k = 1 × 10 −13 m 2 , 1 × 10 −14 m 2 , 1 × 10 −15 m 2 ) so as not to overestimate discharge changes. As a result, coastlines with higher permeability (i.e. coarse-grained unconsolidated sediment) may experience greater discharge changes and permafrost thaw than those projected in this study. More targeted field and modeling studies are needed to assess the impact of local or regional heterogeneity on Arctic CGD.

Conclusions
This study identifies CGD as a likely and increasingly important pathway for the flow of freshwater, and thus carbon, nutrients, and trace metals to the Arctic coastal ocean in a rapidly changing climate. However, results highlight a complex interplay between land-sea hydraulic gradients, warming driven talik formation, and permafrost distribution. This interplay impacts the magnitude and direction of CGD changes, as well as the location and salinity of discharge. Warming-driven talik formation enhances land-sea connectivity and the delivery of freshwater and associated solutes to the coastal zone, yet SLR rate mediates where and how much discharge occurs. Changes in discharge will impact coastal carbon dynamics and biogeochemistry, with the potential to mobilize previously sequestered carbon to the atmosphere or coastal ocean.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.4211/hs.f0f53cbd821842cb817b8c 60c94c4e6f.
is for descriptive purposes only and does not imply endorsement by the U.S. Government.