Upstream ﬂow effects revealed in the EastGRIP ice core using a Monte Monte-Carlo inversion of a two-dimensional ice-ﬂow model

. The Northeast Greenland Ice Stream (NEGIS) is the largest active ice stream on the Greenland Ice Sheet (GrIS) and a crucial contributor to the ice-sheet mass balance. To investigate the ice-stream dynamics and to gain information about the past climate, a deep ice core is drilled in the upstream part of the NEGIS, termed the East Greenland Ice-Core Project (EastGRIP). Upstream ﬂow effects (cid:58)(cid:58)(cid:58) can introduce climatic bias in ice cores through the advection of ice deposited under different conditions further upstream, and are particularly strong at (cid:58) . (cid:58)(cid:58)(cid:58)(cid:58) This (cid:58) is (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) particularly (cid:58)(cid:58)(cid:58)(cid:58) EastGRIP due to its location inside an ice stream on 5 the eastern ﬂank of the GrIS. Understanding and ultimately correcting for such effects requires information on the atmospheric conditions at the time and location of snow deposition. We use a two-dimensional Dansgaard–Johnsen model to simulate ice ﬂow along three approximated ﬂow lines between the summit of the ice sheet (GRIP) and EastGRIP. Isochrones are traced in radio-echo-sounding images along these ﬂow lines and dated with the GRIP and EastGRIP ice core (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) ice-core chronologies. The observed depth-age relationship constrains the Monte Carlo Monte-Carlo method which is used to determine unknown model 10 parameters. We calculate backward-in-time particle trajectories to determine the source location of ice found in the EastGRIP ice core and present estimates of surface elevation and past accumulation rates at the deposition site. Our results indicate that increased snow accumulation with increasing upstream distance is predominantly responsible for the constant annual layer thicknesses observed in the upper part of the ice column at EastGRIP and the inverted model parameters suggest that basal melting and sliding are important factors determining ice ﬂow in the NEGIS. The results of this study form a basis for applying 15 upstream corrections to a variety of ice-core measurements, and the inverted model parameters are useful constraints for more sophisticated modelling approaches in the future. Post-depositional deformation of isochrones observed in radio-echo-sounding (RES) images along ﬂow lines provides information on ice-ﬂow dynamics and can be used to reconstruct past and present ﬂow characteristics. In this study, we use a vertically two-dimensional Dansgaard–Johnsen model to simulate the propagation and deformation of isochrones along three approximated ﬂow lines between the ice-sheet summit (GRIP) and EastGRIP. A Carlo Monte-Carlo method is used to determine the unknown model parameters by minimizing the misﬁt between modelled and observed data. The latter in- cludes the depth of isochrones observed in RES images along the ﬂow lines and a parameter α sur representing the sum of the horizontal strain rates deduced from satellite based surface velocities. of present-day EastGRIP ﬂow

ice temperatures (Salamatin et al., 1998). Processes such as vertical thinning of the ice column and firn densification are also influenced by upstream effects and have consequences on the annual layer thicknesses (Dahl-Jensen et al., 1993;Rasmussen et al., 2006;Svensson et al., 2008) and the age difference between ice and the enclosed air (Herron and Langway, 1980;Alley 55 et al., 1982). Upstream effects in the EastGRIP ice core are expected to be particularly strong due to the fast ice flow in the upstream area (57 ma -1 at EastGRIP, Hvidberg et al., 2020), the strong gradient in the accumulation rate across Greenland's main ice ridge (Burgess et al., 2010), and the increasing elevation towards the central ice divide (Simonsen and Sørensen, 2017). The correction of these effects in the EastGRIP ice core is necessary to interpret the ice core ::::::: ice-core measurements within the climatic context and requires information on the conditions at the time and location of snow deposition. 60 Post-depositional deformation of isochrones observed in radio-echo-sounding (RES) images along flow lines provides information on ice-flow dynamics and can be used to reconstruct past and present flow characteristics. In this study, we use a vertically two-dimensional Dansgaard-Johnsen model to simulate the propagation and deformation of isochrones along three approximated flow lines between the ice-sheet summit (GRIP) and EastGRIP. A Monte Carlo ::::::::::: Monte-Carlo method is used to determine the unknown model parameters by minimizing the misfit between modelled and observed data. The latter in-65 cludes the depth of isochrones observed in RES images along the flow lines and a parameter α sur representing the sum of the horizontal strain rates deduced from satellite based surface velocities. From the modelled velocity field, we calculate particle trajectories backwards in time to infer the source location of ice found in the EastGRIP ice core and estimate the accumulation rate at the time of snow deposition. The source characteristics presented here form a basis to correct for upstream effects in various chemical and physical quantities of the EastGRIP ice core. These corrections are important to remove any climatic bias 70 in ice-core measurements which are currently analyzed and will become available in the future. The inverted model parameters give insight into basal properties and ice-flow dynamics along the flow lines and can be used to constrain more sophisticated numerical models of the NEGIS.

Data and methods
Snow layers deposited at the surface of ice sheets are buried with time and are deformed as a consequence of ice flow. While 75 these isochrones can be observed in RES images today, the ice flow ::::::: ice-flow characteristics which shaped them are generally unknown. This is a typical geophysical inverse problem and can be formulated as d = g(m), where the function g(m) represents the ice flow ::::::: ice-flow model linking the model parameters (m) with the observed data (d). A variety of inverse methods exist to find the model parameters which reproduce the observed data within their uncertainties. Here, we use a Markov Chain Monte Carlo :::::::::: Monte-Carlo : method to determine the unknown parameters of a two-dimensional ice flow :::::: ice-flow : model by 80 minimizing the misfit between modelled and observed isochrones and strain rates. This allows us to reconstruct the ice flow ::::::: ice-flow characteristics in the past and to determine the flow trajectories of the EastGRIP ice.
In the coming sections we describe the data and methods underlying our results according to the work flow illustrated in Fig. 2. In Sect. 2.1 to 2.4 we explain how the isochrone depth-age relationship constraining the Monte Carlo ::::::::::: Monte-Carlo method was obtained. This involves the selection of RES images approximating the EastGRIP flow line (Sect. 2.1), extending 85 the existing chronology of the EastGRIP ice core to the current drill depth (Sect. 2.2), and the tracing and dating of isochrones in the RES data (Sect. 2.3). In Sect. 2.4 the ice flow ::::::: ice-flow model is described in detail and in Sect. 2.5 we elaborate on the Monte Carlo :::::::::: Monte-Carlo method used for parameter sampling. The section numbers are displayed in the corresponding steps in Fig. 2.

EastGRIP flow lines 90
Determining the exact flow line through the EastGRIP ice core ::::::: ice-core site is important to understand the flow history of the survey area and enables us to reconstruct the location where the ice from the ice core was deposited at the ice sheet ::::::: ice-sheet surface. For this, we use high-resolution satellite-based surface velocity products (e.g. Joughin et al., 2018;Gardner et al., 2020;Andersen et al., 2020, see supplementary material Fig. S1) to calculate the upstream flow path. Minor uncertainties and bias in these data products affect along-flow tracing and lead to deviations between flow lines derived from different velocity maps.

95
These deviations become more pronounced with increasing distance from the starting point, as the uncertainties propagate along the line and in general become larger in slow-moving areas of the ice sheet (Hvidberg et al., 2020). Due to the small bias, we consider the DTU_SPACE (Andersen et al., 2020) line the most likely current flow line (Fig. 1b). Yet, there is no evidence that the present-day velocity field was the same in the past. A slight shift in the NEGIS shear margins or the central ice divide, for instance, would have a large effect on the velocity field and, hence, the determination of the flow line of the EastGRIP ice 100 remains ambiguous.
RES data reveal the internal structure of glaciers and ice sheets and provide valuable information on the ice flow ::::::: ice-flow characteristics, particularly when recorded parallel to the ice flow. The electromagnetic waves used in RES are sensitive to contrasts in dielectric properties of the medium they propagate. In ice sheets, these contrasts arise through density variations in the uppermost part of the ice column (Robin et al., 1969), changes in the crystal orientation fabric (Harrison, 1973), and 105 impurity layers such as volcanic deposits (Paren and Robin, 1975). The latter is the most common reflector type below the firn (Millar, 1982;Eisen et al., 2006), and because it is related to layers deposited over a relatively short period of time, most internal reflection horizons (IRHs) detected by RES can be considered isochrones.   Simonsen and Sørensen, 2017;Greene et al., 2017) and the outline of the study area. The NEGIS appears as a distinct feature in the surface velocities (Joughin et al., 2018). It extends from the central ice divide to the northeastern coast, where it splits up into the three marineterminating glaciers 79N Glacier, Zachariae Isbrae and Storstrømmen Glacier. (b) The present-day EastGRIP flow line is derived from the DTU_SPACE surface velocity product (Andersen et al., 2020). Due to the limited availability of radar data along the flow line, we construct three approximate flow lines through a combination of various radar products (profile A-C) between GRIP and EastGRIP. Flow line B and C lack data in the centre of the profiles, marked as a dashed line. The downstream parts of line A and B comprise the same radar profile, which crosses the southern shear margin around 82 km upstream of EastGRIP.  Figure 2. Workflow of the applied steps leading to the results described in Sect 3. The main steps are described in Sect. 2.1 to 2.5 and marked with the corresponding numbers in the figure: The observed data (d obs ) constraining the Monte Carlo :::::::::: Monte-Carlo method consists of the αsur calculated from the ice surface velocities and the isochrone depths along the flow lines. The latter is obtained by approximating the EastGRIP flow line with selected RES images (Sect. 2.1), extending the EastGRIP chronology to the current drill depth (Sect. 2.2) and subsequent tracing and dating of isochrones (Sect. 2.3). The iterative Monte Carlo :::::::::

Initialize model parameters
Monte-Carlo sampling process is illustrated in the grey box (Sect. 2.5) and includes data simulation by a Dansgaard-Johnsen ice flow :::::: ice-flow : model described in Sect. 2.4. connect to two different RES profiles. Line B remains relatively close to the flow direction of the DTU_SPACE line but has a 120 wide data gap in the centre of the profile. In line A, this problem is circumvented by using a radar profile connecting directly to GRIP, which deviates from the observed surface flow field by more than 15 degrees at some locations. Profile C follows the NEGIS trunk all the way to the central ice divide and connects to GRIP over the ice ridge without crossing the shear margin.
Similar to flow line B, flow line C contains a substantial data gap between the onset region of the NEGIS and the central ice divide.

125
To avoid uncertainties related to the proximity of the model boundaries, the flow lines were extended more than 50 km beyond EastGRIP and have a total length of 422 (line A), 421 (line B) and 480 km (line C). To account for any differences in surface elevation or topography between RES data from different years, the ice surface reflections of the radar profiles were aligned to the surface elevation from the Arctic DEM (digital elevation model, Porter et al., 2018). The bed topography in the data gaps of the profiles was derived from the BedMachine v3 data set (Morlighem et al., 2017). The validation of our modelling results and the correct dating of isochrones requires a reliable depth-age scale. The Greenland Ice Core Chronology 2005 (GICC05, Rasmussen et al., 2006;Andersen et al., 2006;Svensson et al., 2006) is based on annual layer counting in various Greenland ice cores. It has been transferred to GRIP and other deep drilling sites in Greenland by synchronizing the ice cores with each other using horizons of e.g. volcanic origin 135 Seierstad et al., 2014). The upper 1,383.84 m of the EastGRIP ice core were drilled between 2015 and 2018, and synchronized with the NorthGRIP ice core in previous work (Mojtabavi et al., 2020).
By 2019, the ice-core drilling progressed down to 2,122.45 m, allowing us to extend the existing time scale from 15 ka to 49.9 ka b2k (thousands of years before 2000 CE). As part of the present study, we identified common isochrones between EastGRIP, NorthGRIP and NEEM to transfer the GICC05 chronology to the part of the EastGRIP record which is not yet 140 synchronized. This involved the same methods applied to NEEM by Rasmussen et al. (2013) and to the upper 1,383.84 m of EastGRIP by Mojtabavi et al. (2020). The isochrones chosen for synchronization purposes are mainly volcanic eruptions, which are registered as brief spikes in the electrical conductivity measurements (ECM, Hammer, 1980). The search of common ECM spikes was performed manually with a strong focus on finding patterns of similarly spaced eruptions rather than single and isolated events. The Matlab program 'Matchmaker' was used to visualize long data stretches and to evaluate the quality 145 of the match ). An iterative multi-observer protocol was applied to reduce problems with confirmation bias and to ensure the reproducibility of the match.
A total of 138 match points were identified between 1,383.84 m and 2,117.77 m, adding to the previously known 381 match points. The match points between EastGRIP and the other two cores are shown in Fig. 3, representing all the volcanic tie points.
The GICC05 chronology was transferred to EastGRIP by linear interpolation of depths between the match points. The age of As in earlier similar work (e.g Rasmussen et al., 2013;Seierstad et al., 2014), very few match points were observed in the stadials, most clearly seen in Fig. 3 in the long stadial stages of GS-2 and GS-3. The sparse volcanic signals within stadial 155 periods should not be attributed to a diminished global volcanic activity but rather to increased deposition of alkaline dust that neutralizes volcanic acid, caused by the prevailing colder and drier climatic conditions (Rasmussen et al., 2013). The largest distance between match points was observed across GS-2 and GS-3 and spans about 162 m of EastGRIP ice.

Tracing and dating of isochrones
The depth-age relationship from ice core :::::: ice-core : chronologies can be extended in the lateral plane by tracing and dating of 'picking tool'. The algorithm is based on calculating the local slope in each pixel of the RES image, and layers are traced automatically between two user-defined points. Starting from each of these points the algorithm walks along the steepest slope 165 towards the other point. Subsequently, the two lines are weighted by distance to their starting point and combined to one layer.
The number of picks required for thorough tracing depends on the data quality and reflector strength.
The total depth uncertainty (z t ) was calculated as where the depth uncertainty introduced during the picking process (z p ) is estimated to be 10 m. The uncertainty related to the 170 radar range resolution (z rr ) of the corresponding RES image is defined as where k is the window widening factor of 1.53, c is the speed of light, B is the radar bandwidth and 3.15 is the dielectric permittivity of ice. time scale. In doing so, local irregularities were smoothed out by averaging the depth over ± 250 m around the trace closest to the ice core ::::::: ice-core location. Because the EastGRIP ice core has not reached the bed yet, we extrapolated the time scale at EastGRIP with two IRHs observed below the current borehole depth to obtain a tentative depth-age relationship between 2,117.77 m and the expected bed depth of 2,668 m.
The total age uncertainty (ã t ) was estimated by following the approach described in MacGregor et al. (2015), where takes into account the age uncertainties associated with the time scale (ã c , equivalent to 0.5 MCE), the radar range resolution (ã rr ), and the layer picking process (ã p ). The uncertainties related to the range resolution are estimated with where a c is the ice-core age from the GICC05 time scale. Similar to Eq. (4),ã p is estimated with The chosen isochrones show distinct patterns which could be identified in all RES images and allowed us to trace isochrones across disruptions and data gaps. Comparison of the isochrone depths at the ice-core locations obtained from different RES images permits to assess ::::::::: assessment :: of the quality of the tracing procedure. The high resolution of the radar images recorded in 2018 facilitates isochrone tracing, and the EastGRIP depths obtained from the two different AWI radar profiles agree : to : within 190 1.5 m. At GRIP, the discrepancy between isochrone depths obtained from three different radar profiles can be up to 30 m, which is slightly above the combined depth uncertainty related to the picking process and the resolution of the RES images. Lower range resolution and signal-to-noise ratio in older RES data introduce bias in isochrone identification, and although distinct isochrones were chosen, a miscorrelation between IRHs recorded by different radar systems can not be entirely excluded.
Moreover, the CReSIS profiles do not precisely intersect at GRIP and deviate from each other. The radar traces closest to GRIP 195 are thus found at slightly different locations for the three RES images, which explains the higher discrepancy of radar layer depths.
The isochrone dating was conducted for each profile individually, and the obtained depths, ages and uncertainties were averaged over the three lines (Table 3). The deepest non-continuous layer which could be identified at EastGRIP is found at a depth of 2,360 ± 11 m and is estimated to be 72,400 ± 1,306 years old. The layer depths of the continuously traced IRHs range 200 from 421 ± 11 m to 2,152 ± 11 m at the EastGRIP location, corresponding to ages of 3,498 ± 94 to 51,920 ± 1,240 years b2k.
Reflectors 1-9 were deposited during the Holocene. The remaining reflectors are found in ice from the Last Glacial Period from which reflector 10 and 11 can be attributed to the onset of the Younger Dryas and the Bølling-Allerød. The relation between the GRIP and EastGRIP depths of the traced IRHs fits well with the GICC05 time scale (Mojtabavi et al., 2020;Rasmussen et al., 2014), and the ages obtained from the two drill sites agree within the uncertainties. We note that the layer dating at

205
EastGRIP consistently leads to younger ages than the dating at GRIP, which is a likely consequence of inaccuracies related to the transformation between ice-core and radar depths.
Due to computational reasons, we did not use all 16 layers for the Monte Carlo :::::::::: Monte-Carlo : inversion but picked ten isochrones with approximately equal vertical spacing, and used the EastGRIP ages for our simulation of layer propagation.
The layers used for the Monte Carlo :::::::::: Monte-Carlo : simulation are indicated in bold in Table 3 and plotted with a consistent 210 color-code in Fig. 4, Fig. 5 and Fig. 7, representing the corresponding ages. Table 3. Characteristics of the traced isochrones connecting the GRIP and EastGRIP ice-core sites. Displayed depths and ages are the average over the three flow lines. Depth uncertainties include the uncertainty related to the picking process and to the radar range resolution. Age uncertainties are related to the GICC05 time-scale uncertainties and isochrone depths. Figure 7d illustrates the depth and climatic context of these layers in the EastGRIP ice core, identified with the corresponding layer numbers. The bold layers and the EastGRIP ages were used for the Monte Carlo :::::::::: Monte-Carlo inversion and are illustrated with a consistent color-code in Fig. 4, Fig. 5  2.4 Ice flow :::::: Ice-flow : model A full simulation of ice flow in the catchment area of the NEGIS is a highly under-determined problem (Keisling et al., 2014), lacking geophysical, climatic and ice-core data, some of which will become available in the future. Simpler models do not solve the problem in detail and are thus computationally much cheaper. Hence, limited but still useful information can be obtained and deformation of internal layers along approximated flow lines between the ice-sheet summit (GRIP) and EastGRIP. The simplicity of the model makes it well suited for the Monte Carlo ::::::::::: Monte-Carlo method due to its few model parameters, the 220 allowance for large time steps, and because it has an analytical solution (Grinsted and Dahl-Jensen, 2002). The model assumes ice incompressibility and a constant vertical strain rate down to the so-called kink height (h) below which the strain rate decreases linearly. Basal sliding and melting are included in the model, and the ice-sheet thickness (H) is assumed to be constant in time.
We consider a coordinate system where the x-axis points along the approximated flow line, the y-axis is horizontal and 225 perpendicular to the flow line, and the z-axis indicates the height above the bed. The horizontal velocities parallel (u ) and perpendicular (u ⊥ ) to the profiles are described by Grinsted and Dahl-Jensen (2002) as: 230 where u ,sur and u ⊥,sur are the surface velocities parallel and perpendicular to the profile, and the basal sliding factor f bed , is the ratio between the ice velocity at the bed and at the surface.
Ice flow in the vicinity of an ice stream is affected by lateral compression and longitudinal extension, in particular across the shear margins of the NEGIS. We thus introduce α = ∂u ∂x + ∂u ⊥ ∂y as the sum of the horizontal strain rates. Due to ice incompressibility we can write α + ∂ω ∂z = 0, where ω symbolizes the vertical velocity. The x and y dependency in Eq. (6 -7) 235 only relates to the surface velocity, such that α sur represents the horizontal dependency in the equations and can be calculated from the ice surface velocities.

240
The boundary conditions for the vertical velocity at the bed (ω bed ) and surface (ω sur ) are whereḃ is the positive basal melt rate,ȧ is the positive accumulation rate, and E bed and E sur are the bed and surface elevations respectively. From Eq. (8) we derive the following expression for the modelled α sur : . (11) Following Grinsted and Dahl-Jensen (2002) and Buchardt and Dahl-Jensen (2007), we adjust the accumulation rates and surface velocities to the climate conditions of the corresponding time with a scaling factor ξ(t): We use the oxygen isotope δ 18 O record from NorthGRIP (Andersen et al., 2004) due to its high temporal resolution, and periods and are defined as (Grinsted and Dahl-Jensen, 2002;Buchardt and Dahl-Jensen, 2007): To simulate the propagation of ice particles deposited at the surface of the GrIS, Eq. (6) and Eq. (8)  Like in most geophysical inverse problems, many different combinations of model parameters can explain the observed data 265 equally well within the range of their uncertainties and therefore, a non-unique solution does not exist. Probabilistic inverse methods consider many different models and describe them in terms of their plausibility, rather than finding one possible solution. This makes it :::: these ::::::: methods : particularly well suited for nonlinear problems, where the probability density in the model space typically shows multiple maxima (Mosegaard and Tarantola, 1995).
Monte Carlo :::::::::: Monte-Carlo : methods are based on a random number generator which allows the sampling according to the 270 target probability distribution in an efficient way. The grey box in Fig. 2 illustrates the iterative sampling process of the Metropolis algorithm (Metropolis et al., 1953) used here: Starting from an initial model (m 0 ), a random walker explores the model space and proposes new models (m new ) which are accepted with a certain probability (P accept ). This way of importance sampling avoids unnecessary evaluation of model parameters in low-probability areas (Mosegaard, 1998).
To estimate the initial accumulation rateȧ 0 , we integrate Eq. (8) (see Appendix A) and obtain the following depth-age where t and z represent the isochrone age and height above the bed, respectively. The accumulation rateȧ is determined with a curve-fitting function, using at least 5 isochrones younger than 10 ka at each point along the flow line. The initial kink height (h 0 ), basal sliding (f bed,0 ) and basal melt rate (ḃ 0 ) are scaled with the normalized surface velocity (û sur ) as follows: where e 1 = 0.4, e 2 = 0.8 and e 3 = 0.03, and the initial value for c w and c c is assumed to be 0.15 and 0.10 respectively.
In each iteration a new model m new is proposed as where m 0 is the initial model and A contains the perturbation amplitude of the corresponding model parameter. The vector q defines the random walk in the multidimensional model space and solely depends on the preceding step. In each iteration, i, one model parameter, j, is randomly selected and perturbed as

290
where r indicates a random number between 0 and 1, and p regulates the maximum step length per iteration of the selected parameter type. To achieve a good performance of the Monte Carlo :::::::::: Monte-Carlo : algorithm, the values of A and p (shown in Table 4) are chosen such that the acceptance ratio for the individual model parameters lies between 25 % and 75 %.
The quality of the proposed model is evaluated by the function S(m) which describes the misfit between the modelled and observed data (see Appendix B). The new model (m new ) is accepted with the acceptance probability (Metropolis et al., 1953) 295 where m old is the last accepted model and the likelihood function is defined as L(m) = e −S(m) .
To ensure that parameter sampling is occurring in a physically reasonable range, the a priori probability distribution is assumed to be uniform within the following intervals: The sampling intervals are based on expected values of the corresponding parameter: The initial accumulation rate obtained from the radar stratigraphy is considered to be quite trustworthy but because the local layer approximation is not justified in the 305 survey area (Waddington et al., 2007) we allow the accumulation rate to deviate by 0.02 ma −1 . The kink height is limited to the ice sheet :::::::: ice-sheet thickness, the basal sliding fraction is allowed to deviate by 30 % from the initial model, and the upper limit of the basal melt rate is based on values suggested at EastGRIP by a recent study (Zeising and Humbert, 2021).
In their initial phase, Markov Chain Monte Carlo :::::::::: Monte-Carlo : methods move from the starting model towards a highprobability area where the target distribution is sampled. To avoid sampling during this so-called burn-in phase, the first 1×10 6 310 accepted models are discarded. Since only one parameter is perturbed at a time, successive models are highly correlated. To obtain a distribution of independent models, only every 1000th accepted model is saved. The sampling is continued for 6×10 6 iterations in total.  The flow-line characteristics and model parameters for each flow line are summarized in Fig. 4. The radar profiles with the observed and modelled isochrones are displayed as a function of the distance from the EastGRIP ice core :::::: ice-core : location.
Particle trajectories were calculated from the simulated velocity field with the mean model parameters and indicate the source location of ice found at the modelled isochrone depth in the EastGRIP ice core. The isochrones and particle trajectories are illustrated with the same color-code as in Fig. 5 and Fig. 7f, indicating the corresponding age. The horizontal strain rates (ε xx , ε yy andε xy ) were obtained from the MEaSUREs Multi-year v1 surface velocity components (Joughin et al., 2018)   IRHs were traced (thin solid lines) in RES images and simulated (dashed lines) with a two-dimensional Dansgaard-Johnsen model (a, e, i).
From the modelled velocity field, we calculated particle trajectories (thick solid lines) backwards in time to obtain estimates of the source location for specific depths in the EastGRIP ice core. The colors of the lines indicate the age of the isochrones and the respective time of snow deposition and are identical to the color code in Fig. 5 and Fig. 7. The horizontal strain rates at the surface were calculated from the MEaSUREs Multi-year v1 (Joughin et al., 2018) surface velocities (b, f, j). The mean and standard deviations of the sampled model parameters accumulation rate, kink height, basal melt rate and basal sliding (c, d, g, h, k, l) were obtained from a Monte Carlo :::::::::: Monte-Carlo inversion by reducing the misfit between observed and simulated data. All panels are aligned at EastGRIP and the x-axis indicates the distance from the borehole location. indicates that the modelled isochrone depth is overestimated which happens in particular for the deepest isochrone towards the end of flow line A and B. As in Fig. 4 and Fig. 7, the color code represents the age of the corresponding isochrone.
sliding factor at EastGRIP, and the accumulation rate in flow lines B and C at EastGRIP. The climate parameter c w is found to be 0.10 ± 0.005 for all flow lines. The obtained value for parameter c c is 0.14 ± 0.003 for flow line A and B, and 0.16 ± 0.004 for flow line C respectively. The histograms of c w and c c can be found in the supplementary material, Fig. S3.

Ice origin and ice-flow history
From the modelled velocity field, we calculate particle trajectories backwards in time (Fig. 4) which give insight into the source location and flow history of ice found at a certain depth in the EastGRIP ice core, and allow us to determine the accumulation rate during its deposition (Fig. 7e). Due to the higher velocities in the ice stream, the ice source location in the upper 1,600 m of the ice core lies further upstream for flow line C compared to flow line A and B. For deeper ice, this trend is reversed, as the 365 velocity along flow line C drops below the velocity of line A and B (Fig. 7a). A similar effect manifests itself in the upstream elevation, where higher velocities along flow line C result in higher elevations in the upper part of the ice column, which is compensated by a flatter topographic profile for ice deeper than 1,400 m (Fig. 7b).
From the model-inferred in situ accumulation rates,ȧ m , and annual layer thicknesses, λ m , we calculate the ice-core thinning function γ: The thinning function increases nearly linearly with depth in the Holocene and shows a considerable decrease in the Younger Dryas and enhanced thinning in the Bølling-Allerød. In the glacial part of the ice core, the thinning function fluctuates be- Table 5. Essential quantities for upstream corrections for selected depths of the EastGRIP ice core. The upstream distance, elevation and past accumulation rates,ȧpast, describe the characteristics of the source location and the conditions during ice deposition.ȧpresent represents the corresponding present-day accumulation rates at the source location. All quantities are averages over the three flow lines and the uncertainties represent the maximum deviation from the mean. tween interstadials and stadials. The shift between the three lines results from the slightly different depth-age relationship and isochrone misfit obtained from the three profiles. We combine the thinning function with the annual layer thicknesses observed 375 in the EastGRIP ice core, λ obs , to estimate past accumulation ratesȧ past : We find that the accumulation rate at the deposition site increases from ∼0.12 ma −1 to a maximum of 0.249 ma −1 for ice at a depth of 912 m, which was deposited approximately 7,800 years b2k. We note that the constant annual layer thicknesses observed in the upper 900 m of the EastGRIP ice core (Mojtabavi et al., 2020) coincides with the spatial pattern of increasing 380 accumulation along the flow line with increasing upstream distance (Fig. 4c,g,k and Fig. 7d,e). Ice between 900 m and 1,400 m is characterized by the transition from the Holocene into the Last Glacial Period with decreased accumulation rates in the Younger Dryas and a peak during the Bølling-Allerød (Fig. 7e). The accumulation rate at the deposition site for older ice varies between 0.02 ma −1 during stadials and 0.196 ma −1 during interstadials. The atmosphere in the glacial period was in general colder and dryer, and hence, accumulation rates were typically lower than today (Cuffey and Clow, 1997). However, due to  Table   3 are displayed with the same color index as in Fig. 4 and Fig. 5 and labeled with the corresponding layer number.
the upstream flow effects, the ice from the interstadials could have been deposited under higher accumulation rates than are observed at the EastGRIP site today.
The variations in the past accumulation-rate between the three flow lines result from both, the varying along-flow accumulation pattern and different upstream distance of the source location. The spread between the three models provides important uncertainty estimates. The average deviation from the mean accumulation rates is 3.9 % in the Holocene and 20 % in the Last

390
Glacial Period. The largest spread between the three flow lines is 68 % observed at a depth of 2,411 m. We remark that, due to missing direct information on the annual layer thicknesses, accumulation rates below the current borehole depth of 2,122.45 m are based on tentative estimates and must be treated accordingly.
The accumulation rates of ∼0.21-0.23 ma −1 at GRIP and ∼0.1-0.13 ma −1 at EastGRIP obtained in this study agree with field observations (Dahl-Jensen et al., 1993;Vallelonga et al., 2014), and the low standard deviations point towards a robust solution. In profiles A and B we observe ∼20 % lower accumulation rates inside the ice stream than outside. This agrees to 405 some extent with Riverman et al. (2019), who found 20 % higher accumulation rates in the shear margins compared to the surrounding, although our observations are not confined to the shear margins only. Regularization on the accumulation rate was necessary in our model to avoid unrealistic strong fluctuations along the flow lines.
The bed topography and bed lubrication have a considerable effect on ice-flow parameters. Flow over bed undulations affect the elevation of internal layers due to variations in the longitudinal stresses within the ice (Hvidberg et al., 1997) and is often 410 reflected in the surface topography (Cuffey and Paterson, 2010). If the bed is 'sticky', i.e. the basal sliding is small, the ice is compressed along the flow direction while vertically extended (Weertman, 1976), and IRHs are pushed upwards. At a slippery bed, the opposite is the case, resulting in along-flow extension of IRHs which leads to thinning and thus decreasing distance between the IRHs. Keisling et al. (2014) argued that major fold trains existing independently of bed undulations can be explained by variations in the basal sliding conditions. This is, for instance, observed across shear margins, where local, steady state 415 folds are formed as a response to the basal conditions (Keisling et al., 2014;Holschuh et al., 2014). In flow line A, we observe similar 'fold-trains' on a larger scale downstream of a substantial bed undulation (100-200 km upstream of EastGRIP) and the resulting high basal melt rates and sliding fraction and the low kink heights drag the layers down in the attempt of matching the observed synclines. These strongly deformed isochrones predominantly appear in parts of the flow lines which deviate from the observed surface velocity direction by more than 15 degrees. We thus argue that they are out-of-the-plane effects and that 420 the isochrones along the ice flow ::::::: ice-flow direction are less strongly deformed. Accordingly, it is questionable that the ice in the EastGRIP ice core has experienced such deformation and that the high local basal melt rates are trustworthy. The fact that these folds are not reproduced very well by the model does therefore not put any constraints on the usefulness of our results for upstream corrections.
The NEGIS differs from other ice streams in Greenland and Antarctica through the lack of clear lateral topographic con-425 straints and high ice-flow velocities reaching exceptionally far inland. The positioning of the shear margins of the NEGIS are most likely strongly interconnected to the subglacial water system and the substrate and morphology of the bed Franke et al., 2021a). The vast amount of ice mass is added to the NEGIS by entering the ice stream through the shear margins (Franke et al., 2021a), resulting in a compressional stress regime perpendicular to the ice stream. The sudden increase in the kink height at around 60 km upstream of EastGRIP pushes the isochrones upwards, similar to the effect of 430 lateral compression. The distribution of available melt-water :::: melt ::::: water and a soft, deformable bed facilitate sliding and thus, ice flow :::::: ice-flow : acceleration at the NEGIS onset . Evidence of a locally enhanced geothermal heat flux and basal ice at the melting point has been presented by e.g. Fahnestock et al. (2001) and MacGregor et al. (2016), and bed lubrication through melt-water production seems to be one of the driving mechanisms for rapid ice flow in the onset region of the NEGIS (Smith- Johnsen et al., 2020b). Our results support these previous findings in the following way: (1)  While it is commonly accepted that the NEGIS is initiated by a locally enhanced geothermal heat flux (e.g. Fahnestock et al., 2001;Alley et al., 2019), the magnitude thereof and the resulting hydrological conditions of the bed are still highly debated.
Previous studies using simple strain-rate models in combination with Holocene radar stratigraphy indicate basal melt rates 445 of 0.1 ma −1 or higher in the vicinity of EastGRIP Keisling et al., 2014;MacGregor et al., 2016).
However, the accuracy of these findings is limited since the local layer approximation (Waddington et al., 2007) is not valid in the surrounding of the NEGIS (Keisling et al., 2014;MacGregor et al., 2016). Remarkably high basal melt rates of 0.16-0.22 ma −1 are also suggested by a recent study (Zeising and Humbert, 2021) using an autonomous phase-sensitive radio-echo sounder (ApRES) at EastGRIP. Melt rates in these order of magnitudes would either require an unusual high geothermal heat 450 flux exceeding the continental background Bons et al., 2021) or an additional heat source (Zeising and Humbert, 2021). Alley et al. (2019) discussed the interactions between the GrIS and the geothermal anomaly, presumably caused by the passage of Greenland over the Iceland hot spot (Lawver and Müller, 1994), and hypothesized that an exceptionally unsteady and inhomogeneous geothermal heat flux underneath northeast Greenland could arise through perturbations of the mantle stress regime caused by ice-sheet fluctuations.

455
Our results indicate basal melt rates at EastGRIP between 0.05 and 0.1 ma −1 but higher values of up to 0.2 ma −1 are obtained further downstream. However, the depth of the oldest modelled isochrone tends to be overestimated in this part of the flow lines (Fig. 5), indicating that the basal melt rate is overestimated. Ice flow :::::: Ice-flow : parameters at a certain location affect the isochrones directly above and further downstream and since EastGRIP is near the end of the radar lines, the information constraining the isochrone depths is limited, leading to overall lower parameter resolution than further upstream.

EastGRIP source location and upstream effects
The source region of ice in the EastGRIP ice core extends over more than 300 km upstream. Holocene ice characterizes the upper 1,244 m of the ice core and has been advected up to 197 km. The climatic conditions during the last 8 kyr remained nearly constant with similar accumulation rates as today. However, due to increasing precipitation towards the central ice divide, ice from the past 8 kyr was deposited under increasingly higher accumulation rates with increasing age (Table 5). Our results

465
indicate that this upstream effect happens to compensate for the vertical layer thinning and results in the constant annual layer thicknesses observed in the upper 900 m of the EastGRIP ice core (Mojtabavi et al., 2020). One possible conclusion of this peculiar observation is that snow depositions must have been advected from far enough upstream to allow the compensation of vertical thinning by increased accumulation rates in the source location. This gives reason to the hypothesis that ice flow ::::::: ice-flow velocities in the past 8 kyr must have been similarly fast as today, and that, therefore, the NEGIS has likely been active 470 during this time. However, we believe that RES images and estimates of present-day accumulation rates along the EastGRIP flow line are necessary to evaluate this hypothesis further.
Between 8 ka b2k and the beginning of the Holocene, accumulation rates decreased at the deposition site due to progressively colder and dryer climatic conditions (Cuffey and Clow, 1997) as we go further back in time and transition into the GS-1. The most recent Glacial Period extends from 119,140 to 11,703 years b2k (Walker et al., 2009) and is characterized by 475 Dansgaard-Oeschger events, abrupt transitions between cold stadial and relatively mild interstadial periods (Dansgaard et al., 1982;Johnsen et al., 1992) causing the oscillations in the accumulation rates. Ice from the Last Glacial Period was deposited between 197 and 332 km upstream from EastGRIP. The basal ice at EastGRIP could be more than 100 ka old which, according to our models, has been deposited within about 50 km from the ice divide under conditions similar to those at NorthGRIP and GRIP.

480
Ice that is entering the NEGIS must somehow penetrate the shear margin, which is an important characteristic of ice flow in ice streams and might have left an imprint on the crystal fabric and texture of ice extracted at EastGRIP. Our modelling results along flow line A and B indicate that ice below 239 m in the EastGRIP ice core has passed the shear margin 82 km from EastGRIP around 1.8 ka b2k. Slightly enhanced annual layer thicknesses observed in the ice core at a depth of 230 m (Fig. 3) seem unrelated to short-term warmer and wetter climate and might thus be an effect of enhanced accumulation across the shear 485 margin, supporting our results.
Our results show surface elevations at the deposition site which are up to 500 m higher than EastGRIP at the corresponding time. Assuming a normal thermal and pressure gradient, this implies that ice was deposited under up to ∼3.25 • C colder temperatures and up to 45 hPa lower pressure than conditions found at the borehole location at the time of deposition.

490
The most relevant limitation of this study arises from lacking radar data parallel to the flow field in the upstream area of EastGRIP. The approximated flow lines deviate from the present-day surface flow field in some parts by more than 15 degrees, which presumably introduces out-of-the-plane effects. Data gaps encumbered isochrone tracing and restricted the Monte Carlo :::::::::: Monte-Carlo : method due to missing information in those areas. The correlated parameters in the Dansgaard-Johnsen model lead to a vast amount of possible solutions, and the fact that the observed data can be reproduced by our model does not 495 prove the validity of the assumed parameters and the physical interpretation thereof. This becomes for instance evident at the GRIP ice core ::::::: ice-core site, where our results indicate basal sliding of up to 30 %, while the drilling project showed that the bed is frozen at the ice sheet ::::::: ice-sheet summit (Dahl-Jensen et al., 1998). The apparant basal sliding might thus represent deformation within a soft bed material rather than actual sliding of the ice over the bed. The spatial and formal resolution of the obtained model parameters is limited, in particular towards the end of the profiles due to limited constraining information 500 further downstream.
By introducing the parameter α, our model accounts for lateral compression and extension on a first degree order, but does not capture the full complexity of the flow field across the shear margins. While these play an essential role in the ice-flow dynamics of the NEGIS  and are likely to have left an imprint on the ice found in the EastGRIP ice core, the full simulation of the flow field is not attempted for the purpose of upstream corrections. This would require more complex, 505 3D numerical ice flow :::::: ice-flow : models which are computationally more expensive and thus not suitable for the Monte Carlo :::::::::: Monte-Carlo : method applied here. Moreover, due to the lack of constraining radar data the information gain in terms of the source characteristics and upstream effects of such a 3D model would be modest.
The elevation of the source location was determined solely from the present-day ice-sheet surface elevation and did not take into account past fluctuations in the ice-sheet thickness. In general, surface elevation changes are relatively minor in the interior 510 areas of central Greenland (Marshall and Cuffey, 2000;Letréguilly et al., 1991). Yet, Vinther et al. (2009) found that the GRIP elevation might have been up to 200 m higher during the early Holocene than today. We did not take into account changes in the ice thickness due to the large uncertainties which would be introduced, particularly in the Last Glacial Period. Our estimates on the surface elevation of the source location must thus not be interpreted as absolute values but rather as relative changes with respect to the surface elevation of the EastGRIP site at the corresponding time.

515
Lacking data and a general understanding of ice-sheet flow far back in time put up additional constraints, and due to the relatively recent discovery of the NEGIS (Fahnestock et al., 1993), little is known about its evolution in the past. Observations of surface elevation and ice-flow velocities imply that the downstream end of the NEGIS has entered a state of dynamic thinning after at least 25 years of stability (Khan et al., 2014). However, it is not clear for how long the NEGIS has been active and how its catchment geometry changed over time.
The assumption of a constant flow field throughout the past 100 kyr is thus the best 520 currently available, but potentially inaccurate, estimate of the past flow regime.
Our results do not give clear evidence on which of the flow lines gives the best results for upstream corrections. Since the present-day EastGRIP flow line is likely located somewhere between flow line A and C, our results can be interpreted as the outer boundaries and we consider the average over the three flow lines the best estimate for the upstream flow characteristics with the corresponding model spread as uncertainties. 525 We traced isochrones in RES images along three approximated EastGRIP flow lines connecting the EastGRIP and GRIP drill sites. A two-dimensional Dansgaard-Johnsen model was used to simulate the propagation of isochrones along these flow lines.
The simplicity of the model allowed us to invert for the ice-flow parameters accumulation rate, basal melt rate, kink height and basal sliding fraction, which give limited but helpful insight into basal properties and ice-flow dynamics and can be used to 530 constrain large-scale ice-sheet models.
On the basis of our modelled two-dimensional velocity field, we calculated particle trajectories backwards in time to determine the deposition site of ice found in the EastGRIP ice core. We present estimates of the upstream distance, surface elevation and accumulation rate at the time and location of ice deposition. This is valuable and necessary information for interpreting ice-core measurements, and to separate past climate variability from non-local imprints introduced by upstream effects. Our 535 studies show that spatially increasing accumulation rates with increasing upstream distance along the flow line are mainly responsible for the constant annual layer thicknesses observed for the last 8 kyr in the EastGRIP ice core.
The lack of radar data along the EastGRIP flow line is the biggest limitation of this study. None of the three simulated flow lines accurately represents the present-day flow field but can be regarded as upper and lower limits framing the upstream effects. The acquisition of further radar data along NEGIS flow lines in the future would thus provide more accurate and valuable 540 insights into the flow history of the EastGRIP ice and the NEGIS.
The RES data recorded by AWI will be available by Jansen et al. (2020) and described by (Franke et al., 2021b). The extended EastGRIP time scale, our derived and approximated flow lines and an extended version of Table 5 will be available on https://www.iceandclimate.nbi.
ku.dk/data/ and in the supplementary material to this paper.
Notation z height above bed x, y direction parallel, perpendicular to the flow profilẽ z t ,z p ,z rr total, picking related, radar related isochrone depth uncertainty k radar window widening factor c speed of light B radar band width a t ,ã p ,ã rr ,ã c total, picking related, radar related, ice core :::::: ice-core related age uncertainty a c ice core :::::: ice-core : age u (z), u ,sur flow-line-parallel velocity at depth, at the surface u ⊥ (z), u ⊥,sur flow-line-perpendicular velocity at depth, at the surfacê u sur normalized surface velocity along the flow line ω(z), ω sur , ω bed vertical velocity at depth, at the surface, at the bed f bed , f bed,0 basal sliding fraction, initial basal sliding fraction h, h 0 kink height, initial kink height H ice thickness α sum of horizontal strain ratesε xx +ε yy E bed bed elevation above sea level E sur surface elevation above sea level a,ȧ m ,ȧ c ,ȧ past ,ȧ present positive accumulation rate, Monte Carlo inferred, ice core :::::::::: Monte-Carlo :::::::: inferred, ::::::: ice-core inferred, past, pres ḃ positive basal melt rate where l runs through the 10 layers, n z and n α are the number of observed isochrone depths and α sur along the flow lines, and M z = d mod,z − d obs,z σ z and M α = d mod,α − d obs,α σ α . (B2)

560
The matrix M z describes the misfit between modelled (d mod,z ) and observed (d obs,z ) isochrones, and the vector M α is the misfit between modelled (d mod,α ) and observed (d obs,α ) α sur . The data uncertainty σ z is the maximum depth uncertainty of 13 m, and the uncertainty on α sur (σ α ) is assumed to be 10 % of the maximum observed α sur . The factor 1000 is a tuning parameter to ensure the acceptance ratio remains between 25 % and 75 %.
Author contributions. DDJ and TAG designed and carried out the study. AG developed the code used for isochrone tracing. AG and CSH derived the EastGRIP flow lines. DJ was co-investigator for the AWI radar survey and acquired the data in the field. SF processed the radar data obtained during the EGRIP-NOR-2018 AWI flight campaign. SOR, GS and TAG synchronized the EastGRIP ice core with NorthGRIP and NEEM and extended the time scale to the current drill depth. TAG prepared the manuscript with the contribution of all co-authors.
Competing interests. The authors declare that they have no conflict of interest.
of Manitoba) and China (Chinese Academy of Sciences and Beijing Normal University). We acknowledge the use of data and data products from CReSIS generated with support from the University of Kansas, NASA Operation IceBridge grant NNX16AH54G, NSF grants ACI-1443054, OPP-1739003, and IIS-1838230, Lilly Endowment Incorporated, and Indiana METACyt Initiative. We also acknowledge the use of the CReSIS toolbox from CReSIS generated with support from the University of Kansas, NASA Operation IceBridge grant NNX16AH54G, and NSF grants ACI-1443054, OPP-1739003, and IIS-1838230. We thank the crew of the research aircraft Polar 6 and system Engineer

580
Lukas Kandora for their work during the AWI flight campaign 2018 and express our gratitude to John Paden and Tobias Binder, who helped with the data acquisition during the AWI survey. Sune Olander Rasmussen and Giulia Sinnl gratefully acknowledge the Carlsberg Foundation for supporting the project ChronoClimate. This research was enabled in part by computing facilities and support provided by WestGrid (www.westgrid.ca) and Compute Canada Calcul Canada (www.computecanada.ca). The scientific colour maps 'hawaii' and 'roma' (Crameri, 2020) are used in this study to prevent visual distortion of the data and exclusion of readers with colourvision deficiencies .