Supplementary Material Radiostratigraphy reflects the present-day, - PDF document

Supplementary Material Radiostratigraphy reflects the present-day, internal ice flow field in the ablation zone of western Greenland Caitlyn Florentine * , Joel Harper, Jesse Johnson, Toby Meierbachtol * Correspondence: Caitlyn Florentine, caitlyn.florentine@umontana.edu 1 Supplementary Methods: Ice Deformation This section outlines details on instantaneous strain, cumulative strain, principal strain, and total strain calculations. 1.1 Ice Velocity and Strain Rate The modelled velocity field is comprised of three components: longitudinal ( 𝑣 ), transverse ( 𝑤 ), and vertical ( 𝑥 ). The model domain is on a Cartesian coordinate system where longitudinal is east-west, transverse is north-south, and vertical is up-down with respect to sea level. Normal strain rates are defined as the spatial gradient of ice flow in the direction of ice flow, so that normal strain rates in the 𝜖𝑣 𝜖𝑤 longitudinal direction are defined as 𝜁̇ 𝑦𝑦 = 𝜖𝑦 , in the transverse direction as 𝜁̇ 𝑧𝑧 = 𝜖𝑧 , and in the 𝜖𝑥 vertical direction as 𝜁̇ 𝑨𝑨 = 𝜖𝑨 . Shear strain rates along a vertical plane extending in the longitudinal 1 𝜖𝑣 𝜖𝑥 direction are defined as 𝜁̇ 𝑦𝑨 = 2 ( 𝜖𝑨 + 𝜖𝑦 ) , along a vertical plane extending in the transverse 1 𝜖𝑤 𝜖𝑥 1 𝜖𝑣 𝜖𝑤 direction as 𝜁̇ 𝑧𝑨 = 2 ( 𝜖𝑨 + 𝜖𝑧 ) , and along a horizontal plane as 𝜁̇ 𝑦𝑧 = 2 ( 𝜖𝑧 + 𝜖𝑦 ) (Hooke, 2005). These were computed on the entire numerical mesh from the Meierbachtol et al. (2015) study, using the finite element solver software FEniCS. The deformation tensor of ice is symmetric (van der Veen, 2013). Thus these six components of normal and shear strain rates provided the full strain rate tensor. 1.2 Finite and Cumulative Strain Approximating strain by taking instantaneous strain rate on a time step (i.e. 𝜁̇ 𝑦𝑧 𝛦𝑢 ) does not accurately resolve large strains (Cuffey and Patterson, 2010), such as those incurred between the ice divide and the margin. Therefore, we calculated finite strain from deposition to emergence. Finite, or logarithmic strain, is defined as (Hooke, 1998): 𝜁 = ln (1 + 𝜁̇ 𝛦𝑢) (I) We compute finite strain on 1 year time steps along P1 (Fig.3), and sum to find cumulative strain.

Supplementary Material 1.3 Principal Cumulative Strain We solved for the eigenvalues and eigenvectors of the cumulative strain tensor to solve for the principal magnitude and direction of cumulative strain at each time step along the modelled particle path. 1.4 A Metric of Total Strain: Natural Octahedral Unit Shear Cumulative normal strain components ( 𝜁 𝑦𝑦 , 𝜁 𝑧𝑧 , 𝜁 𝑨𝑨 ) are used to calculate cumulative natural octahedral unit shear (Hooke, 1998), which is a useful metric of total strain: 1 2 + (𝜁 𝑧𝑧 − 𝜁 𝑨𝑨 ) 2 + (𝜁 𝑨𝑨 − 𝜁 𝑦𝑦 ) 2 ] 2 2 𝛿 𝑝𝑑 = 3 [ (𝜁 𝑦𝑦 − 𝜁 𝑧𝑧 ) (II) 2

Supplementary Material 2 Supplementary Methods: Mass-conserving Ablation 2.1 Relation to Canon The vertically-integrated continuity relation is (Cuffey and Paterson, 2010): 𝜖𝐼 𝜖𝑢 = 𝑏̇ − 𝛼 ∙ 𝑟 ⃑ (III) 𝜖𝐼 Where 𝜖𝑢 represents the rate of ice thickness change, 𝑏̇ represents the surface mass balance, and the ̅𝐼 , where the full ice thickness remaining term represents divergence in ice flux. Ice flux is 𝑟 ⃑ = 𝑉 ̅ = 〈 𝑣 (bed to surface) is H , and the vertical average of ice flow is 𝑉 ̅ 𝑦 , 𝑣 ̅ 𝑧 〉 . Expanded, this flux divergence term is: 𝜖𝑣 ̅ 𝑧 𝜖𝑣 ̅ 𝑦 𝜖𝐼 𝜖𝐼 𝛼 ∙ 𝑟 ⃑ = 𝐼 ( 𝜖𝑦 + 𝜖𝑧 ) + 𝑣 ̅ 𝑦 𝜖𝑦 + 𝑣 ̅ 𝑧 𝜖𝑧 (IV) Plugging this expansion back into the first equation, and assuming ice flow is aligned with the x- coordinate, yields another commonly-used expression of the continuity equation (der Veen, 2009): 𝜖𝑣 ̅ 𝑧 𝜖𝐼 𝜖𝑣 ̅ 𝑦 𝜖𝐼 𝜖𝑢 = 𝑏̇ − [𝐼 ( 𝜖𝑦 + 𝜖𝑧 ) + 𝑣 𝑦 𝜖𝑦 ] (V) Part of the first bracketed term on the right hand side can be simplified using the assumption of 𝜖𝑣 ̅ 𝑧 𝜖𝑣 ̅ 𝑦 𝜖𝑣 ̅ 𝑨 incompressibility, i.e. ( 𝜖𝑦 + 𝜖𝑧 ) = − 𝜖𝑨 = −𝜁̇̅ 𝑨𝑨 . The last bracketed term represents advection of ice thickness gradients. Now we apply this canonical equation to the geometry represented by emergent layer radiostratigraphy in Figure 2. Rather than consider the full ice thickness, we consider a fraction of the ice column, so 𝐼 becomes 𝜃 . We consider surface mass balance averaged along the particle path, ̅ . In the near-surface, where our analysis is restricted, we assume that 𝑉 ̅ = 𝑉 𝑡𝑔𝑑 , thus 𝑏̇ becomes 𝑏̇ and that subsequently −𝜁̇̅ 𝑨𝑨 = −𝜁̇ 𝑨𝑨_𝑡𝑔𝑑 . We consider vertical strain averaged along the particle path, thus −𝜁̇ 𝑨𝑨_𝑡𝑔𝑑 becomes −𝜁̇̅ 𝑨𝑨 , where now the overbar refers to a horizontal average. In a steady state 𝜖𝐼 system, time-dependent rates of change are zero so 𝜖𝑢 = 0 . We consider a steady state because 𝜖𝐼 𝜖𝑢 in the study area are small (Csatho et al., 2014). These substitutions yield: observations of 𝜖𝜃 ̅ − [𝜃(−𝜁̇̅ 𝑨𝑨 ) + 𝑣 𝑦 0 = 𝑏̇ 𝜖𝑦 ] (VI) We consider discrete rather than continuous changes in geometry along the particle path thus 𝜖𝜃 becomes Δ 𝜃 and 𝜖𝑦 becomes Δx. The change in fraction of the ice column is equal to the negative of the initial fraction of the ice column, because 𝜃 1 = 𝜃 at the start point, and 𝜃 2 = 0 at the point of emergence, i.e. Δ 𝜃 = 𝜃 2 − 𝜃 1 = - 𝜃 . We consider the average velocity along the particle path length 𝜖𝜃 −𝜃 𝜃 thus, 𝑣 𝑦 = 𝑤̅ . Assigning this final substitution, the term 𝑣 𝑦 𝜖𝑦 reduces to 𝑤̅ 𝛦𝑦 which is − 𝛦𝑢 . Rearranging terms yields the equation used in our method: 𝜃 ̅ + 𝜁̇̅ 𝑨𝑨 𝜃 − 𝛦𝑢 = 𝑏̇ (VII) 3

Supplementary Material 2.2 Velocity Data Surface ice velocities derived from interferometric synthetic aperture radar (InSAR) satellite data for the winter of 2007 to 2008 are available for the study area at 500 m grid cell resolution (Joughin et al., 2010). We used these observational data to define the average ice speed ( 𝑤̅ ) along each traced radar layer. The direction of ice velocity was parallel to the flight lines. Comparing these velocity data to another satellite-derived velocity dataset (Rignot and Mouginot, 2012), we found the two to differ by <6% along IceBridge flight lines. This is similar to the observational error (~5%; Table 2). We acknowledge that our use of surface velocity observations to define englacial particle travel at depth does not precisely resolve englacial particle speed, but assert that the difference is negligible. Every radar layer analyzed was near-surface, at ice depths <15% of the total ice thickness. Regardless of the mode of motion (i.e. sliding or deformation), ice at such a near-surface depth travels at speeds similar to surface ice (Cuffey and Paterson, 2010). 2.3 Vertical Strain Rate Data A vertical strain rate field is given by the conservation of mass, assuming incompressibility, i.e. 𝜁̇ 𝑨𝑨 = − (𝜁̇ 𝑦𝑦 + 𝜁̇ 𝑧𝑧 ) . When velocity observational errors are propagated through this vertical strain rate calculation, errors are an order of magnitude greater than the vertical strain rate itself, i.e. error ~ 0.01 yr -1 where 𝜁̇ 𝑨𝑨 ~ 0.001 yr -1 . This is a common problem when differentiating discrete data. To address this issue, we smoothed velocity data using a linear convolution at a smoothing window of 4500 m, as it optimized the trade-off between spatial resolution of the velocity field, and magnitude of propagated error through the vertical strain rate calculation (Fig. S2, Supplement). Longitudinal ( 𝜁̇ 𝑦𝑦 ) and transverse ( 𝜁̇̇ 𝑧𝑧 ) strain rates were calculated from the convolved velocity data at a corresponding length scale of 4500 m. 2.4 Ablation Data: K-transect and RACMO In situ melt surface mass balance data are available for 1990-2010 at K-transect sites within the study area (van de Wal et al., 2012). Modelled annual surface mass balances are available at 1 km 2 grid cell resolution from the regional climate model RACMO (Noël et al., 2016). We calculated average ablation rates for RACMO output on a 20 year, modern interval that corresponds with K-transect data (1990-2010). These data are not an exhaustive account or inter-comparison between the many available GrIS regional climate models. However, we assume K-transect data and RACMO output encapsulate ablation measurement uncertainty, spatial variability, and interannual variability characteristic of modern (recent decades) methods and conditions. 2.5 Propagation of Error To honor uncertainties associated with the data used to define each term ( η , Δx , 𝑤̅ , and 𝜁̇̅ 𝑨𝑨 ) in Eq. (1), we propagate observational error according to a Taylor series expansion, given by: 2 𝜖𝐺 𝑓𝑠𝑠𝑝𝑠 = ∑ 𝐺 [ 𝜖𝛽 𝑗 𝛽 𝑗𝑓𝑠𝑠𝑝𝑠 ] (VIII) 𝑗 where F is the function, α is each input variable ( i ) in the function, and error is denoted by subscript. 4