Mobilization of subsurface carbon pools driven by permafrost thaw and reactivation of groundwater flow: a virtual experiment

Permafrost thaw leads to an increase in groundwater circulation and potential mobilization of organic carbon sequestered in deep Arctic sediments (e.g. 3–25 m below surface). Upon thaw, a portion of this carbon may be transported along new groundwater flow paths to surface waters or be microbially transformed or immobilized by in-situ biogeochemical reactions. The fate of thaw-mobilized carbon impacts surface water productivity and global climate. We developed a numerical model to investigate the effects of subsurface warming, permafrost thaw, and resultant increased groundwater flow on the mobilization and reactive transport of dissolved organic carbon (DOC). Synthetic simulations demonstrate that mobilization and groundwater-borne DOC export are determined by subsurface thermo-chemical conditions that control the interplay of DOC production (organic matter degradation), mineralization, and sorption. Results suggest that peak carbon mobilization from these depths precedes complete permafrost loss, occurring within two centuries of thaw initiation with the development of supra-permafrost groundwater flow systems. Additionally, this study highlights the lack of field data needed to constrain these new models and apply them in real-word site-specific applications, specifically the amount and spatial variability of organic carbon in deep sediments and data to constrain DOC production rates for groundwater systems in degrading permafrost. Modeling results point to key biogeochemical parameters related to organic matter and carbon bioavailability to be measured in the field to bridge the gap between models and observations. This study provides a foundation for further developing a physics-based modeling framework to incorporate the influence of groundwater flow and permafrost thaw on permafrost DOC dynamics and export, which is imperative for advancing understanding and prediction of carbon release and terrestrial-aquatic carbon exchange in warming Artic landscapes in the coming decades.


Introduction
Permafrost degradation facilitates greater partitioning of liquid water into the subsurface and enhances groundwater circulation (Bense et al 2009, Young et al 2020, Yao et al 2021, leading to increased groundwater discharge to Arctic rivers (Smith et al 2007, Walvoord and Striegl 2007, St. Jacques and Sauchyn 2009. With permafrost thaw, vast stores of organic carbon sequestered in permafrost are susceptible to mobilization and transformation to greenhouse gases (Schuur et al 2008, Hugelius et al 2014, Olid et al 2022. While past research has mostly focused on vertical carbon release to the atmosphere, thaw-mobilized carbon can also be laterally transported by groundwater flow to surface waterways , O'Donnell et al 2012, Vonk et al 2015, Sun et al 2021. Hence, permafrost loss and subsequent increases in groundwater discharge alter the magnitude of both total carbon release and dissolved organic carbon (DOC) export from terrestrial to aquatic environments (McGuire et al 2010, Connolly et al 2020, Sterte et al 2022. Organic matter reservoirs in permafrost environments generally consist of the vegetation canopy, leaf litter and shallow soils (∼0-3 m) where seasonal biogeochemical cycling of 'young' carbon occurs. Deeper stores of 'old' carbon also exist in sediments containing buried organic matter that can be several tens of meters thick and tens of thousands of years old (Strauss et al 2013, Schuur et al 2015. These 'old' carbon pools are generally not part of the active carbon cycle in a stable climate . However, with climate warming and permafrost thaw, carbon can be reintroduced into the critical zone and active carbon cycle (McGuire et al 2010, Wild et al 2019, as evidenced by observed increases in old carbon contributions from deeper thawed subsurface layers to streamflow (Connolly et al 2020, Olid et al 2020). Global estimates of deep organic carbon (3-25 m below surface) stored in permafrost areas, such as the Yedoma region and Arctic river deltas, range from 210 to 465 Pg (1 Pg = 1 × 10 12 kg) (Schuur et al 2015), which is up to ∼60% of the carbon mass presently stored in the atmosphere (Houghton 2007).
The thawing of organic matter increases temperature-dependent decomposition and production of DOC (O'Donnell et al 2012), and increased groundwater circulation mobilizes DOC (Wickland et al 2018, Connolly et al 2020. However, deeper thaw and longer groundwater flow paths increase groundwater residence time and opportunity for DOC sorption onto unfrozen mineral surfaces (Kawahigashi et al 2006), but higher temperatures also increase microbially-driven mineralization (Striegl et al 2005, Lyon andDestouni 2010). Consequently, the net impact of thaw on DOC mobilization and export from thawing permafrost will be determined by the interplay between these competing mass transport and attenuation processes (O'Donnell et al 2012, Wickland et al 2018. Such mechanisms have been proposed to explain differences in carbon release and DOC export across watersheds with different states of thaw; however, the dearth of field data and limited modeling capabilities for considering cryohydrogeological and biogeochemical processes have limited conclusive assessments of the conditions under which each process dominates. Thus, improved understanding of the effects of groundwater dynamics on carbon mobilization and processes governing carbon export in thawing landscapes is crucial for predicting changes to Arctic carbon cycling .
Herein we develop physics-based numerical simulations that consider coupled groundwater flow, thermal energy transfer with phase change, and reactive solute transport to elucidate time scales and hydrogeological conditions that control permafrost DOC mobilization by groundwater under climate warming. Simulations investigated the production, reactive transport, and export of leached DOC in an evolving groundwater flow system typical of thawing Arctic hillslopes. Simulated DOC dynamics were used to evaluate the relative importance of competing reactive transport processes and expected timescales relevant for permafrost carbon mobilization and export in natural systems.

Methods
The finite element modeling environment FlexPDE (PDE Solutions 2020, www. pdesolutions.com) is used to solve the coupled partial differential equations describing groundwater flow, reactive solute transport, and thermal energy transfer with freezethaw. FlexPDE has been successfully benchmarked and frequently applied for simulating groundwater flow in thawing groundwater systems (e.g. Bense et al 2012, Grenier et al 2018); these models have recently been expanded to investigate solute transport processes in northern groundwater systems (Guimond et al 2021. We further expand on these efforts by incorporating temperature-dependent, reactive solute transport processes specific to carbon mobilization from thawing permafrost.

Water and energy transport
Assuming saturated subsurface conditions with constant density, groundwater flow is described by (Bear 1972): where, h (m) is the hydraulic head, ρ w (kg m −3 ) is groundwater density, g (m s −2 ) is acceleration due to gravity, µ (kg m −1 s −1 ) is dynamic viscosity of water, k (m 2 ) is the aquifer intrinsic permeability, S s (m −1 ) is the porous medium specific storage, and k rw (−) is a relative permeability exponential function that reduces permeability, by up to six orders of magnitude, as pore ice is formed (McKenzie et al 2007). Energy transport and subsurface temperature (T ( • C)) distribution, including the effects of thermal conduction, advection and latent heat transfer, are described by (e.g. Bense et al 2009): where κ e (W m −1 K −1 ) and C a (J m −3 K −1 ) are the effective thermal conductivity and volumetric heat capacity of the sediment/fluid/ice mixture, respectively, q (m s −1 ) is the groundwater (Darcy) flux, and ρ i (kg m −3 ) is the ice density. C w (J m −3 K −1 ) and L i (J kg −1 ) are the volumetric heat capacity and massbased latent heat of fusion of water, respectively. The reduction in unfrozen water-content (θ w (m 3 m −3 )) with freezing temperatures as ice forms is described using an exponential soil freezing curve (McKenzie et al 2007, Devoie et al 2022. Expressions for all constitutive relations and parameter ranges used in equations (1) and (2) for various model scenarios are listed in supplementary table S1.

DOC reactive transport
The leaching of organic carbon into pore water and transformation to DOC involves temperaturedependent production, transformation and transfer of organic carbon that act collectively to add or remove DOC in groundwater before it discharges to surface receptors (Kawahigashi et al 2006). If organic matter stabilized in permafrost is the source of carbon, upon thaw, organic carbon can be partitioned into: (a) immediate release of CO 2 and CH 4 gases, (b) particulate organic carbon, and (c) porewater DOC, the species simulated here. Once leached into groundwater, DOC undergoes microbially-driven mineralization and adsorption/desorption to mineral surfaces (Fan et al 2010). The advection-dispersion equation is a suitable DOC transport model provided that DOC production, mineralization, and sorption processes are considered (Yurova et al 2008). To make our model numerically tractable when coupled to water flow and heat transfer, the following assumptions are employed: (a) a constant zero-order source term represents DOC production (Yurova et al 2008), on the assumption that organic matter (OM) is not limited in the region of the carbon reservoir over the time-scale simulated; (b) a linear equilibrium sorption approach describes the net sorption of DOC to mineral surfaces (Dusek et al 2019); and (c) firstorder kinetics represent mineralization of aqueous and sorbed phases of DOC (Fan et al 2010). In total, DOC processes include advective-dispersive transport with OM degradation to DOC (production), DOC sorption, and DOC mineralization of both aqueous and sorbed phases: where C is the DOC concentration (g m −3 ), D (m 2 s −1 ) is the hydrodynamic dispersion tensor (Bear 1972), v i (m s −1 ) is the linear groundwater velocity (v i = q/θ w ), P (g m −3 s −1 ) is the temperature-dependent zero-order DOC microbial decomposition (production) rate, ρ b (kg m −3 ) is the aquifer bulk density, n (m 3 m −3 ) is the aquifer porosity, S C (g kg −1 ) is the adsorbed concentration of DOC onto the aquifer matrix, and λ a and λ s (s −1 ) are temperature-dependent first-order mineralization coefficients for the aqueous and sorbed DOC phases, respectively. The temperature dependency of DOC production and mineralization terms are represented using a modified Van't Hoff (i.e. Q 10 ) equation (table S2) describing the exponential increase in the rate of biological processes with temperature (Yurova et al 2008). Expressions for all constitutive relations and values for all parameters used in equation (3) are listed in table S2.

Model simulations
Simulations represent archetypal two-dimensional (2D) hillslope configurations including a carbon reservoir as sediments containing a pool of initially unreactive organic matter encased in permafrost. Upon thaw, carbon from this reservoir leaches into groundwater as DOC (figure 1). Similar hillslopescale archetypal simulations have been performed without reactive transport of carbon, to examine the interacting effects of permafrost thaw and groundwater flow ( Hydraulic head along the surface is fixed (h = surface elevation) to represent the position of the water table (figures 1(a) and (b)) as only the saturated zone is considered. Based on the imposed water table slope (figure 1(b)), recharge occurs at the top of the hillslope, while discharge is focused at the lower gradient side of the land, representing direct discharge to a stream plus surrounding wetlands and riparian areas. The model base and vertical sides are closed for fluid flow (figure 1(b)). For heat transport, a 65 mW m −2 geothermal heat flux is applied along the base, typical of the subarctic (Blackwell and Richards 2004), while the vertical sides are closed for heat exchange (figure 1(c)). Temperatures along the top boundary were based on thaw scenarios that considered a mean annual land surface temperature subject to a linear warming trend described later. Seasonal freezethaw was not considered due to this study's focus on deeper carbon and DOC as well as the impacts of decadal to centennial thaw on permafrost extent (Bense et al 2012).
The solute reservoir that exists (figure 1(d)) at the initial condition (t = 0) represents sediments containing a pool of previously stable organic matter from which DOC can be produced upon thaw (Neff et al 2006. Boundary conditions for solute transport included a flux boundary condition along the surface such that any recharging water has a DOC concentration of 0 g m −3 . Along the discharge portion of the hillslope, which evolves over the simulation period, advective DOC fluxes exit the model domain with discharging groundwater, representing DOC export to surface water. The model base and sides are closed to solute mass fluxes (figure 1(d)). Geological, geochemical, and geothermal complexities were intentionally minimized to examine fundamental aspects of groundwater flow, heat transport, and reactive solute transport processes controlling groundwater discharge and DOC mobilization.
Initial conditions for groundwater flow, heat transport, and thus permafrost conditions (figure 1) represented an equilibrium condition (simulated for 1000 years) for a mean annual temperature of −2 • C applied to the domain surface to generate permafrost. The carbon reservoir (figure 1(d)) is present at depths 3-23 m below the ground surface, representing a 20 m zone in the model that assumes the mineral sediment contains significant fractions of organic matter (Dutta et al 2006, Strauss et al 2013 that is initially frozen. Following thaw, i.e. when temperatures in the region rise above 0 • C, this becomes the source of DOC in the system due to organic matter degradation and leaching of DOC. This conservative depth range of the carbon reservoir (3-23 m) was chosen as these depths are where most deep carbon pools in the Arctic are thought to be hosted within unconsolidated sediments since the last glacial maximum (Hugelius et al 2014). From initial conditions, simulations were run for 600 years to allow for complete permafrost thaw. A warming trend of +0.05 • C yr −1 was imposed on the initial mean annual surface temperature for 100 years, followed by a 500 year period at the new stable mean annual temperature of 3 • C. This trend is compatible with the coupled model intercomparison project 5 (CMIP5) representative concentration pathway 5 (RCP8.5) warming projections for Arctic regions (Collins et al 2013) and consistent with current Arctic warming trends (Landrum andHolland 2020, Rantanen et al 2022). Aquifer permeability was highest at the surface and decreased logarithmically with depth (table S1) (Quinton et al 2008). These landscape parameters were based on generalized subsurface conditions for regions in North America and Eurasia hosting large sedimentary deposits of organic-rich silty sediments, e.g. the Yedoma region (Schirrmeister et al 2011). The permeability distribution and water table boundary condition are very conservative from a water flow perspective as they generate annual groundwater recharge rates ranging from 0 to 20 mm for fully frozen to fully thawed conditions respectively. These rates are well below precipitation in high-latitude systems (Serreze andHurst 2000, Sheffield et al 2006) and in line with recharge estimates for the Arctic ((Döll and Fiedler 2008), so as not to impose unrealistically large aquifer flushing rates. A total of 27 simulations were performed with DOC reactive processes (i.e. production, mineralization, and sorption) varying from high to low representative parameter values (table 1).

Choice of DOC parameters
DOC is primarily released into environments through decomposition of plant residues and release from shallow soil carbon pools (Fan et al 2010). Thus, shallow DOC is driven primarily by modern plantderived DOC as near-surface organic carbon undergoes production-leaching-accumulation cycles. However, this study focuses on DOC production and mobilization from deeper stores of 'old' carbon buried in landscapes since the last glacial maximum (Schuur et al 2008(Schuur et al , 2015.
Most studies that have simulated DOC leaching and export have focused on shallow active carbon pools in temperate and boreal forest soils or northern peatlands, which generally have high organic matter content compared to deep carbon pools (Schuur et al 2015). For example, estimates range between 1% and 10% organic carbon by weight in deep permafrost sediments in the Yedoma region, and these tend to be ice-rich where large volumes of ground-ice decreases the organic carbon content by 22% and 50% compared to shallow soil carbon in Arctic permafrost (Dutta et al 2006, Strauss et al 2013, Schuur et al 2015. Furthermore, shallow soil environments are often highly oxygenated, unsaturated (in the vadose zone), warmer, and have much higher primary productivity than deep permafrost carbon, which impacts the structure, lability, and reactivity of the carbon. Thus, organic matter concentrations as well as degradation and DOC production rates in these warmer environments will generally be much higher than permafrost carbon in saturated groundwater systems, which will likely have a large portion of anoxic conditions with concomitant lower DOC production rates.  Weston andJoye 2005, Yurova 2008). Examples of DOC production in permafrost groundwater systems could not be found in the literature, but the large difference measured in marine sediment porewaters (Weston and Joye 2005) suggest that DOC production will be orders of magnitude lower than shallow soil systems. Yurova et al (2008) found DOC production in a permafrost-affected peatland complex to be around 1.16 × 10 −7 g m −3 s −1 . As this was for a permafrost environment, but still oxygenated, this value was used as the high end of the range used (tables 1, S3 and S4). Microbially-driven mineralization rates are less uncertain as measurements have been made in groundwater systems, and saturationdependent mineralization relationships have been developed (Yurova et al 2008, Mei et al 2012. The values used here are for saturated systems. Sorption is controlled by aquifer bulk density and mineralogy, and thus we employed a range of sorption coefficients for silty sediments. Table 1 summarizes all the model runs and associated reactive transport parameters used, with further explanation provided in supplementary tables S3 and S4.

Effects of thaw-driven groundwater flow on timescales of peak DOC mobilization
Results show that permafrost thaw increased groundwater discharge (figures 2(a) and (b) DOC mobilization, consisting of both microbially-driven DOC mineralization and groundwater-driven DOC export, begins 60 years after the onset of warming (figure 2(c)), or 20 years after permafrost loss commenced ( figure 2(b)). Due to the assumption of continued organic matter degradation and thus DOC production, and uncertainty in the amount of deep carbon, only the first 200 years of the results were used to analyze DOC mobilization on an order of magnitude basis. Figures 2(d)-(g) show the range of mobilization fluxes (mineralization + export) and mass attenuation (stored DOC) results based on the different reactive parameters employed (table 1). The subsurface DOC mass in both the aqueous and sorbed phases ranged from 10 1 to 10 5 g m −1 due to the assumption of equilibrium sorption (figures 2(d) and (e)), with DOC mass in the sorbed phase slightly higher due to its lower rate of mineralization (Kalbitz et al 2005).
Rates of total DOC mobilization from the onset of thaw over the first 200 years ranged from 10 −6 to 10 −4 g s −1 m −1 (figure 2(c)). Of this, the majority of the mobilized DOC is via mineralization (figure 2(f)), as rates of permafrost DOC export over the first 200 years were smaller but exhibited a much larger flux range of 10 −11 -10 −5 g s −1 m −1 (figure 2(e)). Figure 3 shows the cumulative carbon mobilized over 200 years, which ranged from 5 × 10 2 to 3 × 10 6 g m −1 for mineralized DOC and 5 × 10 −2 -3 × 10 4 g m −1 for groundwater-facilitated exported DOC. While DOC export is generally much lower than DOC mineralization, in many scenarios (e.g. C7, C13 and C16 in figure 3) DOC export is not insignificant and is as high as DOC mineralization for other scenarios (e.g. C1, C10 and C19 in figure 3).
This large range of mobilization and export resulted from only a three order-of-magnitude variation of DOC reaction parameters (production, mineralization, and sorption), which were used as proxies for varying aquifer biogeochemical characteristics controlling how leached permafrost DOC is transformed and transported through the subsurface (section 2.4). Despite this variability, the timeline of export was consistent across simulations, i.e. within the first 200 years of simulation, DOC mobilization rates reached within 90% of their simulated maximum (figure 2(c)), well before complete permafrost loss (blue line in figure 2(b)). Figure 4 shows the simulated hydraulic heads, icesaturation, DOC concentrations (normalized to each scenario's maximum concentration), and temperatures at (a) 100 years, (b) 200 years, and (c) 500 years. The figure illustrates how hydrogeological and permafrost thaw responses controlled the evolution of the subsurface flow and fate of thawmobilized DOC, including: (a) top-down thaw into the organic matter reservoir producing and mobilizing DOC, and (b) development of a coupled groundwater flow, thermal regime, and reactive transport system controlling DOC fate. As surface temperature increases and permafrost thaws, pore-ice melts into liquid water which increases the hydraulic conductivity of the thawed region by up to 6 orders of Figure 3. Range of cumulative carbon mobilization (export + mineralization) on a log scale from time 0 to 200 years for each simulation scenario C1 to C27 (listed in table 1) along with their respective production rate (top), mineralization rate (middle) and sorption coefficients (bottom), with the text color indicating high (red), medium (green), and low (blue) values. Scenarios with values of DOC export lower than 1 g were considered negligible and not plotted. * Vertical scale is slightly higher than other panels. magnitude and creates a gradually expanding permeable supra-permafrost groundwater flow system (see Lamontagne-Hallé et al 2020 and references within). The increase in temperatures above 0 • C allows microbially-driven degradation of previously frozen organic matter and DOC production, while thaw simultaneously creates permeable subsurface flow paths that allow the release, transport, and reaction of DOC in the system as well as export to surface water. These coupled water, heat and solute transport processes collectively drive the release of previously sequestered carbon into the active hydrological and carbon cycle. During the warming period (first 100 years), groundwater discharge increases with the deepening permafrost table and development of the supra-permafrost aquifer ( figure 4(a)). Top-down thaw exposes formerly frozen sediments to above-zero temperatures, making the organic matter within them biogeochemically available for DOC production and reactions (figures 4(a) and (b)).

Coupled subsurface transport processes driving DOC mobilization
Groundwater circulation is initially limited to this upper portion of the aquifer (figures 4(a) and (b)). At approximately 500 years, permafrost is thawed to the point that deeper flowpaths begin to develop (figure 4(c)) and within 40 years after this time all permafrost has disappeared ( figure 2(b)) and the full hillslope-scale flow system develops.
In all simulations, DOC mobilization rates reached 90% of their maximum value within the first 200 years (figure 2(c)). However, by 200 years, the groundwater discharge has only reached 50% of its maximum ( figure 2(b)). While these results would depend on the initial permafrost thickness, this demonstrates that simulated peak DOC fluxes were driven more by rising temperatures and newly created permeable flow paths in the supra-permafrost aquifer than the increase in groundwater discharge from the development of a regional flow system which occurs at a much later stage. The simulated increase in subsurface temperatures show that warming in the DOC reservoir in the last 300 years was minimal compared to the first 200 years which drove the initial activation of DOC production (OM degradation). This also demonstrates that DOC mineralization and export were largely mediated by the thermal regime in the warming supra-permafrost groundwater system ( figure 4(b)), not the activation of deep flow paths ( figure 4(c)). This suggests that the largest export of old carbon from the depths simulated here could occur before full permafrost degradation and reactivation of regional groundwater flow systems. This implies: (a) mobilization of permafrost DOC begins well before substantial permafrost loss, as is being observed currently; and (b) that even if surface temperatures stabilize at the current levels in the near future, substantial increases in DOC mobilization and groundwater-derived export to Arctic streams of old, previously inactive carbon are likely to still occur.
It should be noted that this study employed a 2D archetypical model. The 2D approach made simulations computationally feasible and enabled processbased assessment of the coupled DOC-permafrostgroundwater flow system in a changing climate. However, 2D approaches overlook three-dimensional (3D) water, energy, and solute transport as well as multi-dimensional variability in permafrost extent, solute dispersion, topography, organic matter distributions, and subsurface heterogeneity. However, this investigation provides a first-of-its-kind numerical analysis of coupled groundwater, heat, and solute transport processes driving permafrost carbon release and DOC mobilization and serves as a catalyst for future investigations and model development. It also calls for integration of field data and 3D simulations to capture the multidimensional nature of hydrogeological processes (e.g. Evans et al 2015, Langford et al 2020).

Future work and data needed to link models to observations
The simulation results elucidate how permafrost thaw and hydrogeological reactivation could result in the mobilization of deep pools of carbon from DOC mineralization and export to surface waters. The amount and rates of DOC mobilization depend on organic matter quantity and quality, and site-specific hydrogeological, thermal and biogeochemical aquifer conditions. The range of DOC parameters employed here were designed to represent a plausible range of conditions that drive DOC mobilization across diverse Arctic landscapes and groundwater systems. The range in simulation results (figures 2 and 3) show that the magnitude of permafrost DOC mineralization and export to streams, as well as the ratio of the two fluxes is strongly controlled by subsurface characteristics that mediate temperature-dependent Figure 5. Data needed to bridge gaps between observations and models, and to inform parameterization of the new type of model presented in this study. Blue panels = hydraulic data and parameters, red panels = temperature data and parameters, and green panels = reactive transport data and parameters. Panels with bolded text highlight key data parameters needed for DOC simulations and thicker arrows indicate key linkages between data and model parameters. production, mineralization, and sorption of subsurface DOC. However, there are critical gaps in observational data that limit the applicability of the new type of model employed here to real-world, sitespecific conditions. Figure 5 illustrates the hydrogeological, thermal, and biogeochemical data needed to constrain the model's flow and transport parameters. Of these, the data most lacking are the biogeochemical observations that drive DOC production and mineralization rates (Pi et al 2021, Wright et al 2022. Hydrogeological parameters such as aquifer permeability, initial permafrost thickness ( figure 5 are ground and subsurface temperatures as the thermal regime of permafrost and supra-permafrost groundwater systems drives the depth and timing of permafrost thaw, chemical reaction kinetics, and soil microbial activity. While the variability in hydrological and thermal parameters was not fully explored here, these parameters are more easily constrained with currently available datasets (e.g. Zipper et al 2018, Evans et al 2020, and are not the limiting factors in allowing this model to be applied to real-world settings. Furthermore, the objective of this research was not to evaluate model sensitivities to these parameters (see de Bruin et al 2021 andHamm andFrampton 2021 for examples), but to gain insight into the information on the subsurface biogeochemical conditions needed to improve predictions of the potential response of sequestered carbon to permafrost thaw (figure 5).
In addition to hydrogeological and thermal properties, subsurface biogeochemical data are needed to inform site-specific models and improve predictions of potential responses to permafrost thaw. DOC production is driven by the decomposition of organic carbon present in permafrost sediments, thus measurements of organic matter concentrations below 3 m are essential, yet globally sparse (Hugelius et al 2014). Thus, the amount and spatial variability of deep organic carbon concentrations are among the most important data needed to constrain the response of deep permafrost carbon to climate warming. In this regard, information on the impact of ice-rich permafrost on carbon concentrations is also needed. Just as important as organic carbon content and distribution with depth is the bioavailability of permafrost organic carbon and DOC (e.g. molecular structure and chemical recalcitrance) for which there are even scarcer measurements (Pi et al 2021). DOC mineralization is driven by biogeochemical conditions such as aquifer/soil mineralogy, pH, dissolved oxygen content, redox conditions, and microbial activity. Anoxic versus oxic conditions are likely the most important parameters as they would determine whether DOC is mineralized to CO 2 or CH 4 (Prokushkin et al 2009, Wagner andLiebner 2009). In-situ dissolved gas measurements in permafrost groundwater systems would also constrain mineralization rates used in the type of model presented here.
Previous modeling studies on DOC reactive transport in shallow soil systems (e.g. Yurova et al 2008, Mei et al 2012, Dusek et al 2019 have constrained their production and mineralization rates by calibrating models to DOC concentrations (for production) and greenhouse gas efflux measurements at the ground surface (for mineralization). However, simulation of deep carbon pools as done in this study is more complicated as mobilized carbon from deep stores will mix and interact with shallow mobilized carbon, and thus distinguishing the fluxes from these respective DOC sources is less straightforward. While carbon isotopes can be used to assess old carbon contributions to DOC export (e.g. Gonzalez Moguel et al 2021), temperature-controlled incubation experiments on undisturbed sediment cores taken from permafrost sediments would potentially allow for direct measurements of DOC production and mineralization rates. Incorporating these types of measurements into the type of model presented here would help constrain DOC production, mineralization and export in site-specific applications. The current lack of these data results in considerable uncertainty in the biogeochemical consequences of thawing permafrost . Collection and measurements of these data shown in figure 5 will be useful for characterizing the potential or 'risk' for old carbon release from permafrost landscapes, as well as the fate of such carbon, based on subsurface conditions rather than solely the amount of carbon present.
It is also crucial to note that the simulations presented here do not consider the fate of DOC once mineralized or exported with groundwater discharge, as we cannot distinguish between the specific mineralization pathways or the mineralized DOC interactions with shallow soils, vegetation, and aquatic systems with this modeling framework. Models that accommodate instream and near-surface soil DOC processes are far more common than the subsurface model presented here. Efforts are needed to link subsurface reactive transport models with models that incorporate instream and soil DOC processes to consider the interaction of surface, shallow and deep carbon pools.

Summary and conclusions
This study developed a physics-based model that considered the effects of thaw-activated Arctic groundwater flow systems on the mobilization and reactive transport of deep permafrost carbon. Results show that DOC from deep sediments (>3 m) could be a significant source of carbon to the active Arctic carbon cycle as the climate warms. Simulations further suggest this could occur well before complete permafrost loss and development of deep groundwater flow paths. This implies that a maximum mobilization and export of deep stores of permafrost carbon could occur under current warming rates within the next 2 centuries.
As the simulations presented here only demonstrate first-order controls on reactive transport behavior, considerable uncertainty remains in the range of total DOC mobilization and export given the complexities of biogeochemistry in permafrost settings. This uncertainty provides the impetus for more fieldbased measurements to inform future simulations. However, as shown here, by analyzing the influence of subsurface thermal regimes on biogeochemical activity, important aquifer characteristics can be used to provide assessments of deep carbon vulnerability. Consideration of these subsurface processes will be needed to constrain the Arctic response to permafrost thaw and quantify terrestrial feedbacks to climate change.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.