Horizontal Ice Flow Impacts the Firn Structure of Greenland's Percolation Zone

. One-dimensional simulations of firn evolution neglect horizontal advection from ice flow, which transports the firn column across climate gradients as it is buried by accumulation. Using a suite of model runs, we demonstrate the impacts of horizontal 15 advection on the development of firn density, temperature, and the stratigraphy of melt features through the Greenland ice sheet percolation zone. The simulations isolate processes in synthetic runs, and investigate four specific transects and an ice core site. Relative to one-dimensional simulations, the horizontal advection process tends to increase the pore close-off depth, reduce the heat content, and decrease the frequency of melt 20 features with depth by emplacing firn sourced from higher locations under increasingly warm and melt-affected surface conditions. Preservation of the advected pore space and cold content is strongly dependent upon the depth of meltwater infiltration. Horizontal ice flow interacts with topography, climate gradients, and meltwater infiltration to influence the evolution of the firn column structure; the interaction between these variables 25 modulates the impact of horizontal advection on firn at locations around Greenland. Pore close-off and firn temperature are mainly impacted in the lowermost 20-30 km of the percolation zone, which may be relevant to migration of the lower percolation zone. the percolation zone, however, the stratigraphy of melt features can have an advection derived that should be 30


INTRODUCTION
Summer melting of bare ice, epitomized by stream networks and moulins, represents a relatively small portion of the Greenland Ice Sheet (GrIS) periphery since about 90% of the 35 ice sheet's area is perennially snow covered accumulation zone (e.g., Ettema et al., 2009). A large fraction of the snow covered region also experiences melt ( Figure 1): between 50-80% of the area melted during summers of the period 1958-2009 (Fettweis et al., 2011), for example. Further, the inland extent and duration of melting have demonstrated increasing trends and have frequently established new records (Mote, 2007;Tedesco, 40 2007; Tedesco et al., 2013). Melting of the accumulation zone (i.e., the percolation zone) is therefore an increasingly important aspect of the ice sheet, and so too are the glaciological processes governing the snow/firn interactions with surface climate.
Meltwater from the lower accumulation zone may run off from its point of origin (e.g., 45 Machguth et al., 2016), while at higher elevations the water may simply infiltrate into cold snow and firn to fill underlying pore space, forming ice when it refreezes (e.g., Braithwaite et al., 1994;Harper et al., 2012) or remaining liquid if it does not (e.g., Forster et al., 2014;Humphrey et al., 2012). The capacity of the firn column to accommodate meltwater is dependent on its thermal state, its density, and structure of ice layers. However, many 50 aspects of the processes governing firn's structural and thermal evolution, and whether meltwater is retained, remain unclear. While current model fidelity prevents confident constraint on the amount of melt retained in the percolation zone, existing estimates are that 40-50% of the meltwater generated never escapes (van Angelen et al., 2013;Janssens & Huybrechts, 2000;Reijmer et al., 2012). 55 The percolation zone is a region with relatively high horizontal motion compared to submergence rate (cf. divide regions) (Figure 1). Ice sheet flow displaces the firn column to lower elevation, where it is buried by subsequent winter layers experiencing higher intensity summer melt. Thus, the deep firn column's structural makeup and thermal state 60 results from a climate that varies in both time and space. Ice motion potentially impacts the structure of the firn column including the amount of deep pore space that could absorb meltwater and act as a source of associated latent heat. Further, it may have implications for the interpretation of melt feature stratigraphy within ice cores collected from these regions. 65 Here we investigate the role that horizontal motion plays in driving the structural evolution of the firn layer in Greenland's percolation zone. We utilize previous approaches for modeling firn densification and meltwater infiltration, but we extend our analysis to include horizontal advection of the domain due to ice flow. Climate and ice flow are highly 70 variable around the ice sheet, but sparse observational data, computational limitations, and questions of fidelity surrounding models for meltwater infiltration in firn, prevent us from simulating the entire ice sheet with a high level of confidence. Our purpose here is therefore to test for the importance of ice flow to firn structure in Greenland's percolation zone, because ice flow has largely been overlooked and is not currently included in regional 75 climate model simulations of firn evolution. To explore this process, we focus our investigation on synthetic modeling of isolated processes, four differing transects of the GrIS percolation zone, and partitioning the signal of climate change from an advection signal within an ice core from the percolation zone.

Model Description
Air temperature, accumulation rate, and melt/refreezing processes drive the evolution of the density and thermal structure of the percolation zone's firn column (e.g., Herron and Langway, 1980;Reeh et al., 2005). The spatial gradients in these parameters, coupled with 85 the speed at which the ice moves through the gradients, also impacts evolution of the firn column. We simulate these processes in a thermo-mechanically coupled framework for firn densification and heat transfer, including meltwater penetration and refreezing, that also incorporates horizontal displacement of the firn column due to ice motion. The rate of change of firn densification is described by its material derivative: where ρ is firn density and u is the velocity vector. Many existing densification models are based on the empirical assumption that the proportional density change is linearly related to the change in overlying load (Robin, 1958), and ours is no different. As rate parameters, 95 we use the empirical constants developed by Herron and Langway (1980) for dry snow densification, based upon the relatively simplistic formulation with few tuning parameters and favorable comparison with other densification schemes (Lundin et al., 2017). The total densification rate is therefore: where the critical density ρ 550 kg m , and temperature and accumulation-dependent 100 constants c and c are defined as: (3) In Equation 3, R is the gas constant (8.314 J K -1 mol -1 ), accumulation rate b is in ice equivalent units, and T is the absolute temperature.

105
Firn temperature is simulated using the standard heat transfer equation, including latent heat additions from meltwater refreezing (see Section 2.1.3): where c is the specific heat capacity of ice, taken to be 2100 J kg -1 K -1 , and S reflects latent heat release from refreezing. In Equation 4, K is the thermal conductivity, which we prescribe to be density-dependent following Arthern and Wingham (1998) (K=2.1(ρ/ρi) 2 ).

Melt Infiltration Schemes
Modeling complex and heterogeneous meltwater infiltration in firn with high fidelity remains an outstanding problem of critical importance and is beyond the scope of this 115 project. Our approach is to implement three existing infiltration schemes which vary in complexity and reflect a range of approximations. The first model considers only shallow infiltration, assuming that all meltwater refreezes in the top annual layer (Reeh et al., 2005). The second implements a standard bucket method (Kuipers Munneke et al., 2014;Ligtenberg et al., 2018), allowing meltwater infiltration as far as permitted by thresholds 120 for cold content and irreducible water content (the latter of which is defined following Coléou and Lesaffre (1998)). Meltwater percolates until reaching either a firn layer in which the available liquid water fails to exceed the layer's irreducible water content or the pore close off density; any remaining meltwater runs off instantaneously. The third infiltration model implements a continuum approach (Meyer and Hewitt, 2017), simulating 125 the physics of water flow based on Darcy's Law, and treating both saturated and unsaturated conditions.

Numerical Methods
We implement horizontal advection in modeling exercises using an approach that is as stress enhancement induced by horizontal deviatoric stresses (Alley and Bentley, 1988), but are likely to be negligible around the GrIS percolation zone (Supplemental Material 145 1.1).
In addition, our modeling approach neglects horizontal heat diffusion. To assess the consequences of this, we tested our Lagrangian approach against an explicit 2D model for densification and heat transport including horizontal diffusion in an Eulerian framework.

150
Results from these tests yielded negligibly different results (Supplemental Material S1.2; Figure S1). We therefore continue with our Lagrangian approach, which streamlines computational efficiency and permits flexible implementation of various melt infiltration schemes.

165
The influence of horizontal advection on firn structure at depth is dependent on ice flow speed and spatial gradients in climate forcings (temperature, melt, and accumulation). We conducted an initial test of model sensitivity to each of these variables to understand, in isolation, the influence of changes in these processes on firn structure. We then applied the model to four flowline transects across GrIS' percolation zone representing a spectrum of

Greenland Transects
Our modeling approach was implemented at four test transects ( Figure 1) around the GrIS percolation zone: 1) the well-studied EGIG transect in western GrIS, 2) a transect feeding Jakobshavn Isbrae, 3) the K-transect in southwest GrIS, and 4) a transect extending into

195
Helheim Glacier. These four study profiles were selected to capture a wide variety of ice sheet conditions (Table 1). The inland extent of the percolation zone was selected to correspond with the location where melt, as a fraction of accumulation rate, is equal to the value at Crawford Point ( Figure 1); a choice based on temperature measurements indicating melt infiltration at Crawford Point is insufficient to warm the firn above the 200 annual average temperature . The initial condition for the interior boundary of each transect was the steady state profile based on the values given in Table 1.
Surface velocities along study transects were defined from satellite velocity data (Joughin et al., 2010), and RACMO2.3p2  was used to select 1980-2016 average 205 climate variables ( Figure S3). This time period captures the increase in GrIS melt since the late 20th century (Fettweis et al., 2011). In addition to the transect simulations incorporating horizontal ice flow, in each transect we also completed 1D simulations at 600-1700 locations spaced at annual displacements (calculated from surface velocities) along the flowline. The latter were used for baseline comparisons of the effects of including 210 or neglecting horizontal advection of the firn column.

Impact on Core Stratigraphy
A commonly used metric for quantifying changing climate conditions from firn cores is the annual increment of surface melt, or Melt Feature Percent (MFP) (Graeter et al., 2018; 215 Kameda et al., 1995;Koerner, 1977;Trusel et al., 2018). To investigate the role that horizontal advection can play in MFP records from the percolation zone, we simulated the conditions leading to Crawford Point (69.877°N, 47.0102°W, 1997 m elev), located along the EGIG line. This site is relatively high in the percolation zone; in recent decades the average summer at this site experiences about 15 days of melt (Mote, 2007). Multiple 220 shallow cores have been collected for density and temperature measurements Humphrey et al., 2012), and in 2007 and a deep core was collected and interpreted within the context of GrIS melt history (Higgins, 2012). The site therefore offers an opportunity to assess the role of horizontal advection on interpretation of melt history in a core profiling the full firn column of the percolation zone. 225 We modeled the 2D firn evolution along a flow line beginning 22 km inland, and ending at Crawford Point using datasets for the modern state. This transect extent is chosen to ensure that the simulated conditions at Crawford Point contain no remnants of the initial condition. Ice surface geometry (Morlighem et al., 2017(Morlighem et al., ) and, mean (1980(Morlighem et al., -2016 melt and 230 snowfall values from RACMO2.3p2  were used to determine spatial climate gradients. As with the other Greenland transects (Section 2.3.2), horizontal advection during burial is represented by present-day velocity datasets (Joughin et al., 2010). We employ the (Reeh et al., 2005) model for infiltration to be consistent with past MFP observational studies which assume all annual melt is confined to the corresponding 235 annual layer (e.g., Graeter et al., 2018;Kameda et al., 1995;Trusel et al., 2018). We assume the spatial gradients in input datasets have not changed over the century time scale required to simulate firn conditions at the bottom of the firn column. The validity of this assumption is unknown and perhaps tenuous; our intention, however, is a demonstration of the horizontal advection process constrained by ice sheet conditions. Furthermore, if 240 there are in fact large time changes in gradients, this only increases complexity in the signal created by horizontal advection.

245
Including 2D horizontal advection in simulations of the percolation zone yields greater air content in the firn column and increased depth to pore close off compared to 1D results Adding horizontal advection to simulations also decreases the firn temperature; the temperature profile and temperature at pore close off reflect advected firn from higher, colder conditions. Heat content is strongly influenced by the choice of meltwater routing 270 scheme: for example, under very high accumulation and melt, the bucket method yields deep penetration of water and warmer firn temperature at depth (cf. the 1D case) ( Figure   S2). Steeper topography yields larger along-flow gradients in melt, temperature, and accumulation, causing greater disparities between 2D-advection and 1D-profile simulations. The ice flow speed has potential to impact simulations with 2D-advection, but 275 importantly, the magnitude is strongly modulated by the values and gradients in other variables. Interestingly, while meltwater infiltration and refreezing amplifies the difference between 1D and 2D-advection simulations of firn density, the release of latent heat eliminates advected cold content, thereby reducing the thermal disparities between 1D and 2D simulations ( Figure S2).

Transects
The most significant differences between the 1D and 2D model simulations are along the lowermost 20-30 km of our four sample transects. By including ice flow in these firn simulations, the density decreases by >50 kg m -3 for the EGIG, Jakobshavn, and Helheim 285 transects ( Figure 3; Figure S4; Figure S5), resulting in increases of pore close off depth of up to 8 m, 13 m, and 19 m, respectively. The commensurate impacts on total air content in the firn column can also be large: for example, along the EGIG transect it changes by ~50% in the lower 10 km, and by 7%-25% along the next 10-20 km.

290
The different melt infiltration schemes yield variable impacts. As in the sensitivity testing, the largest impact is with the Reeh et al. (2005) scheme, under which the inclusion of horizontal advection in simulations increases the firn column air content by up to several meters from a 1D simulation (Figure 4). In contrast, the air content changes are almost always less than one meter in the continuum scheme, which allows for deep meltwater 295 penetration.
Local changes in surface slope along the transects both enhance and diminish the impacts of horizontal advection on the underlying firn structure, complicating the 2D geometry of the firn layer within the percolation zone. This is exemplified in both EGIG and Jakobshavn 300 transects ( Figure 3). A flattening of the surface profile near 20 km (EGIG) and 30 km (Jakobshavn) means that snow accumulates under similar climate conditions over many km. This generates a deep firn density structure that is more similar to 1D profiles as the firn is buried and horizontally advected through the transect. These flat regions terminate abruptly with a sharp increase in slope, and therefore climate gradients as well. This 305 results in deep firn with densities that are reduced by 20-30 kg m -3 compared to 1D simulations. In contrast to the EGIG and Jakobshavn transects, the changes to density structure throughout the K-transect are comparatively small. Surface speeds are consistently ~18 -50% of the values along EGIG and Jakobshavn counterparts (Table 1); such slow velocity all but eliminates the impact of ice flow (Figure 3d).

310
The process of horizontal advection generates colder firn temperature profiles. Along the EGIG transect, horizontal advection decreases firn temperatures at the depth to pore close off by 1.0°-1.5° C in the lower 15 km, and by 0.8°-1.0° C in the next 20 km. With the high speeds, steep topography, and heavy melt of the lowermost reaches of Jakobshavn and 315 Helheim transects, firn temperatures were reduced by as much as 3° C by including horizontal advection.

Impact on Core Stratigraphy
At Crawford Point, firn deposited along the first ~5km above Crawford Point reflects a 320 region with very low slope and essentially no climate gradient caused by elevation ( Figure   5A). Consequently, MFP values are relatively constant in the upper ~50m of the firn column, indicating that horizontal advection is inconsequential to firn structure over depths that are equivalent to recent decades at this site ( Figure 5B).

325
Below ~50m depth, however, there is an abrupt inflection to continuously decreasing MFP to the bottom of the core ( Figure 5B). Horizontal advection generates a decline in MFP of approximately 7 percentage points from 50 m depth to the end of the core. Presented in terms of simulated age (rather than depth), this amounts to an apparent reduction in melting of 0.04 percentage points per year that arises from horizontal ice flow and not time 330 changes in climate. As discussed below, this is a non-trivial magnitude when scaled against changes in melt arising from warming climate.

335
The choice of meltwater infiltration scheme has a large effect on the simulated impacts of horizontal advection of firn in the percolation zone and is a key uncertainty in the fidelity of model results. In reality, water moves vertically as a wetting front propagating downward from the surface (Colbeck, 1975), but also by complex and inhomogeneous infiltration processes (Marsh and Woo, 1984;Pfeffer and Humphrey, 1996), and it can be routed 340 horizontally along impermeable ice layers (e.g., Machguth et al., 2016). With so little known about deep infiltration, none of our schemes are likely to be entirely accurate: the Reeh et al. (2005) scheme only allows melt penetration within the annual snow increment which is known to be incorrect, especially low in the percolation zone where melt rates are high (e.g., Humphrey et al., 2012); the continuum model (Meyer and Hewitt, 2017) uses the 345 most complex physics, but has large uncertainties for coefficients of permeability and grain sizes; and, the bucket model (Kuipers Munneke et al., 2014;Ligtenberg et al., 2018) disregards the complex physics governing flow of water through the firn matrix, simplifying the problem to just density and cold content and assuming the flow of meltwater is instantaneous.

350
With horizontal advection of firn tending to move open pore space underneath an increasingly melting surface, the depth/quantity of infiltration is key: the deeper melt penetrates, the more the deep pore space is 'overprinted' by surface melt and not preserved. Alternatively, infiltration that is limited to shallow depths enhances the 355 disparity between deep firn and that nearer to the surface. Our suite of model runs show that, in the lower percolation zone, the choice of infiltration scheme has nearly equivalent impact on the total air content as the incorporation of ice flow.

360
Simulation results demonstrate the changing influence of horizontal advection on firn structure in response to changing climate gradients and ice speed as the firn parcel traverses the percolation zone. Transect modeling indicates that along any given flowline the influence of advection tends to increase towards the lower percolation zone; an intuitive result considering that surface speed and slope (a proxy for climate gradients) 365 both increase substantially relative to the upper percolation zone, and the surface experiences heavy melt (Supplemental Figure S3). Sensitivity testing showed that each of these factors amplifies differences between 2D and 1D representations of deep firn structure.

370
Transect modeling also reveals that ice flow introduces variability in firn structure around the ice sheet, not just along individual flow lines. Surface speeds within the GrIS' percolation zone shown in Figure 1 vary from nearly 0 to more than 1000 m a -1 . The K-Transect, EGIG, and Jakobshavn transects demonstrate the differences in firn structure that can develop regionally across Western Greenland as a result of ice flow patterns (Figure 3).

375
However, the EGIG and Jakobshavn simulations show that advection can also influence firn structure at a more local scale. Despite being separated by just ~40 km, differences in surface speed develop between the two transects in the lower 20-30 km as ice in the Jakobshavn transect accelerates towards the margin ( Figure S3). This results in a simulated firn column that is 5-10 m thicker compared to the nearby EGIG profile at the same 380 elevation. While the local gradients in ice speed are perhaps greater here than nearly anywhere else on the ice sheet, the local and regional differences in our simulated transects illustrate that differences in deep firn structure are likely to exist in regions of the GrIS percolation zone with otherwise similar climate conditions, purely as a result of differences in ice flow patterns and topography that dictate spatial gradients in climate. While the MFP trend at Crawford Point is impacted by horizontal advection, simulated 420 profiles of firn density and temperature have little influence from the advection process relative to lower on the transect (Figure 3). This result may seem counterintuitive.
However, the density and temperature fields evolve over a time-space continuum, whereas the MFP record represents a small time-trend in the occurrence of discrete events.
Occasional thin ice layers, which are easily identified in cores, can have little impact on the Certainly some locations in the percolation zone may yield ice cores with MFP trends that are not significantly impacted by ice flow. But considering the potential for ice flow to obscure climate trends, a simple procedure for quantifying this effect has utility. If the 435 present ice sheet state (speed, accumulation, and melt rates) is assumed to be constant in time, an apparent climate signal at any core site can be quantified from spatially extensive datasets of the above variables. At a dated core depth corresponding to a time t years before present, the firn parcel originated at a location (x(t)) upglacier from the core location, where x(t) is the integral of the spatially varying velocity ( along the flowline 440 over t years: . The MFP at time (t) can therefore be determined from the accumulation and melt conditions at this upglacier location: .
Equations 5 and 6 can thus be combined to generate a time series of MFP at a core site that is a record of spatially varying climate advected by ice flow; the component that should not 445 be incorrectly interpreted as time-changing climate.

CONCLUSIONS
Elevated horizontal ice flow in the percolation zone compared to ice divides results in a firn column that is not always well represented by 1D models for time-evolving density and 450 temperature. shortcomings of a 1D simulation. The horizontal advection process thus has greatest influence on firn evolution in the lower accumulation zone (e.g., 20-30 km); a nexus of conditions that are likely migrating upward as climate warms, but are also subject to the greatest uncertainty regarding melt infiltration processes.

460
The 2D evolution of firn in the percolation zone is influenced by topography: horizontally invariant firn is generated in flat regions, whereas local hills/swales enhance the 2D influences from horizontal advection. The deeper meltwater penetrates, the more pore space is filled by surface melt and the advected deep pore space and cold content is not preserved.

ACKNOWLDEGEMENTS
Funded by NSF grants 1717241 (Harper, Meierbachtol) and 1717939 (Humphrey), and a Montana NASA Space Grant Fellowship (Leone). The authors thank J. Johnson for discussions, and B. Vandecrux and an anonymous reviewer, whose careful reviews improved the manuscript. This paper has no data to declare. All model simulations are available at https://github.com/um-qssi/Firn_advection. (a) delineated by black contour lines. Thick black lines through percolation zone show study transects, where E is EGIG, J is Jakobshavn, K is K-transect, H is Helheim (see Table  1).

Figure 2.
Example sensitivity test. Modeled differences between 1D and 2D for ice speed 625 using dry firn model (black), Reeh model (red), tipping bucket model (blue), and continuum model (green). Base accumulation is 0.5 m a -1 , approximately the average value of the EGIG transect shown in Figure 1. Each simulation is run with surface slopes of 0.3° (dotted), 0.6° (solid), and 0.8° (dash-dotted) to represent different climate gradients that may exist around the GrIS.