Firn Evolution at Camp Century, Greenland: 1966–2100

Camp Century is an American military base built in 1959 under the surface of the Greenland ice sheet and decommissioned in 1967. Here, we use outputs from RACMO2.3p2 and CanESM2 climate models, adjusted to meteorological observations, and a firn model to simulate the firn density and temperature at Camp Century between 1966 and 2100. The model output is evaluated against an extensive set of firn 3observations and three Representative Concentration Pathways (RCP2.6, 4.5 and 8.5) are considered as future scenarios. Our model suggests that the upper horizon of the Camp Century debris field – observed at a depth of 32 m in 2017 – will continue to be buried by persistent net accumulation over the next eighty years under all RCP scenarios. This horizon depth will be between 58 and 64 m in 2100, depending on the RCP scenario. We estimate a maximum meltwater percolation depth of 1.1 m under all RCP scenarios. We therefore find it extremely unlikely that surface meltwater interacts with the subsurface debris field at Camp Century before 2100 under all RCP scenarios. Camp Century’s future is representative of the firn area in northwestern Greenland, bound to shift from dry snow to a percolation regime. Our model suggests that 10 m firn temperatures at Camp Century will increase from −24.0°C in 1966 to −21.3, −20.0 and −18.6°C in 2100 under the RCP2.6, 4.5 and 8.5 scenarios, respectively. We reveal a previously unknown warm bias in air temperatures simulated at Camp Century by both RACMO2.3p2 and CanESM2 climate models which needs to be accounted for when using these models to predict melt, firn evolution and sea-level contribution of the Greenland ice sheet. We also present novel in situ measurements of firn compaction rates, which indicate that about 25% of firn compaction of the top 62 m of firn occurs below 20 m depth. This highlights the importance of deep-firn compaction measurements for model evaluation and correction of altimetry products.


INTRODUCTION
Camp Century is located in the northwest corner of Greenland's high-elevation firn plateau, at an elevation of 1886 m, approximately 200 km east of Thule Air base (Figure 1).Camp Century was constructed in 1959 by the United States Army Corps of Engineers as a cut-and-cover trench network within the relatively soft and porous firn (Clark, 1965;Mellor, 1969).The base operated year-round until 1964, and then seasonally for an additional three summers, before being permanently abandoned with minimal decommissioning in 1967.Persistent snowfall has now buried the base entirely, leaving nothing visible at the ice-sheet surface.Under the downward advection associated with each year's snow accumulation, the tunnel network has descended through the firn.An ice-penetrating radar survey of Camp Century showed that from an initial depth of 8 m in 1959 the tunnel network was between 45 m and 55 m depth in 2017 (Karlsson et al., 2019).The data suggest that 95% of the subsurface debris field is now located at depths greater than 32 m.This subsurface debris field is approximately circular with a radius of less than 1 km.
Camp Century represents a crucial site for two reasons.First of all, surface mass balance at the site has been projected to shift from net accumulation to net ablation by 2100 under a high-end greenhouse gas emission scenario (Fettweis et al., 2013;Colgan et al., 2016).Consequently, the Government of Greenland raised concerns regarding the potential remobilization of contaminants at Camp Century via meltwater percolation.As a response, the Government of Denmark established the Camp Century Climate Monitoring Program in 2017, led by the Geological Survey of Denmark and Greenland (GEUS), with a goal "to regularly update annual likelihoods of meltwater interacting with abandoned materials at the Camp Century site over the next century" (Colgan et al., 2017).Secondly, Camp Century is located in the accumulation area of the Greenland ice sheet where surface melt has recently increased in frequency and magnitude due to anthropogenic climate change (van den Broeke et al., 2016;Fettweis et al., 2017).In these regions, the firn can refreeze and retain the meltwater generated at the surface each summer and buffer the ice sheet's mass loss and sea-level contribution (Benson, 1962;Pfeffer et al., 1991, Harper et al., 2012).The Greenland firn responds to increased meltwater through either a deeper percolation and refreezing of meltwater (Humphrey et al., 2012), the formation of impermeable near-surface ice layers, with subsequent increased local meltwater runoff (de la Peña et al., 2015;Machguth et al., 2010;MacFerrin et al., 2019), or, in the case where the excess meltwater neither refreezes or runs off, the development of perennial firn aquifers where meltwater remains liquid (Forster et al., 2014;Miège et al., 2016).Due to a long history of meteorological and glaciological observations, Camp Century represents a reference site to document the past, current and future evolution of the firn layer in northwestern Greenland.This study therefore aims at describing the evolution of the firn at Camp Century from 1966 to 2100.
For this purpose, three domains need to be described: 1) the near-surface atmosphere, that can be characterized by standard meteorological variables; 2) the snow surface where energy and mass exchanges take place; and 3) the firn itself.All three domains interact and adapt to each other continuously.The near-surface atmospheric conditions can be monitored by Automatic Weather Stations (AWS) or simulated by climate models.While AWSs give the best estimation of the local weather, they are limited in time, they are subject to frequent malfunctions, and cannot be used for future climate projections.Climate models do not suffer gaps and can be run using scenarios for future climate but are known to deviate from observed weather conditions due to their coarse resolution and inherent simplifications.We consequently adjust four outputs from two climate models to the observations from two AWSs enabling us to describe the present and future meteorology at Camp Century.We use two AWSs covering discontinuously the 1996-2019 period: the GITS station from the Greenland Climate Network (GC-Net, Steffen et al., 1996) and CEN station from the Camp Century Climate Monitoring Program (Colgan et al., 2018).
For climate model outputs, we first use, over the 1966-2019 period, the regional climate model RACMO2.3p2(Noël et al., 2019) forced by re-analyzed climate data, which is based on meteorological observations.These climate data represent our best knowledge of the surface conditions at Camp Century when AWS data are not available.Secondly, the output of an Earth System Model (ESM), CanESM2 (Chylek et al., 2011), is used over the 1966-2100 period.ESMs simulate the entire Earth's climate evolution and, under certain assumptions about planetary evolution, can estimate future climate.As future scenarios, CanESM2 use the Representative Concentration Pathways (RCPs) employed by the United Nation's Intergovernmental Panel on Climate Change (IPCC) in their Fifth Assessment Report (van Vuuren et al., 2011;Myhre et al., 2013).We here use the RCP2.6,RCP4.5 and RCP8.5 scenarios.They project a radiative forcing, enhanced by greenhouse effect, rising to +2.6, +4.5 and +8.5 W/m 2 by 2100, respectively, and associated global mean temperature changes of +1.5, +2.9 and +4.6 °C in 2081-2100 relative to the pre-industrial period (Rogelj et al., 2012).Studies of sea-level budget in each scenario estimated that RCP2.6, 4.5 and 8.5 respectively lead to a global mean sea-level rise of 0.4, 0.47, and 0.63 m by 2100 (Church et al., 2013).The atmosphere interacts with the firn through mass exchange such as snowfall, sublimation or deposition and through the surface energy budget.The surface energy budget describes all the energy fluxes occurring, through the surface, between the atmosphere and the firn: the radiation received and emitted by the surface, the sensible heat flux due to the temperature difference between the air and the surface, the latent heat flux due to the phase change of water vapor at the surface and lastly the energy transferred from the firn to the surface.Various calculation schemes exist to convert meteorological information into surface energy fluxes (e.g., Box and Steffen, 2001;Radić et al., 2017).We here use the surface energy budget model from van As et al. (2005) as used in Vandecrux et al. (2018), Vandecrux et al. (2020a).
Given the surface mass and energy fluxes, the firn model simulates the evolution of firn density and temperature through time.This is done by discretizing the firn into many vertical layers and solving equations for thermomechanical continuity as mass and energy transfer between these layers.Several different firn models are in use around the world today, differing slightly in their formulations of physical firn processes, or in their numerical methods (see reviews from e.g., Steger et al., 2017or Vandecrux et al., 2020b).Most of these models include the meltwater infiltration, refreezing and runoff, the diffusion of heat, the release of latent heat during meltwater refreezing, and the densification of snow into ice.We here use the GEUS multilayer firn model, which was presented and evaluated at various sites in Vandecrux et al. (2018), Vandecrux et al. (2020a), Vandecrux et al. (2020b).
This study therefore evaluates and adjusts four climate model outputs (one deemed most reliable for the historical period and three providing predictions until 2100) to AWS data, calculates for each of them the surface mass and energy budget and uses them as input to a firn model.The firn model outputs for the four forcing datasets is evaluated against a suite of firn measurements and used to describe the recent evolution of the firn from 1966 to 2019.Eventually, the firn model forced by the three future climate predictions is used to describe the potential evolution of the firn to 2100 and assess whether the military waste will interact with meltwater in the next decades.

Historical Review: Firn and Weather Observations at Camp Century
Unique archives of historical firn density and temperature measurements have been collected at Camp Century.The earliest historical measurements date to a 1952 SPIRE (United States Snow, Permafrost and Ice Research Establishment) traverse expedition with joint reconnaissance and scientific mandates that passed the site where Camp Century would eventually be located (Benson, 1960).While this glaciological traverse was repeated annually until 1955, only the 1954 data is readily available today.This dataset includes measurements of firn density and temperature to depths of 8.8 m at "Station 2-20" (Benson, 1962).The next known survey is from a subsequent 1956 SIPRE polar glaciology course.This expedition also measured firn density and temperature at several sites along the over-snow trail originating at Thule, including "Mile 120", but the original data is not readily available today (Roch, 1956;Ragle, 1958).
The United States Army Corps of Engineers (USACE) began drilling a deep ice core at Camp Century in 1964.After two unsuccessful attempts, drilling on the ultimate deep core was started in 1964 and finished in 1966 (Ueda and Garfield, 1968).The near-surface density profile of this core record therefore pertains to 1964, while the borehole temperature profile was measured after drilling in 1966.While original ice temperature measurements are available, these measurements were made at 50-m down-borehole depth increments (Weertman, 1968;Gundestrup et al., 1993).This limits the utility of the 1966 ice temperature profile in examining near-surface firn.While Camp Century was active, the USACE conducted intensive research to characterize firn density and temperature-dependent strength (McCoy and Waterhouse, 1960;Mellor, 1969).The vast majority of these engineering studies, however, applied to the disturbed firn of gallery walls and are therefore not useful for this study.
In 1969, a joint USACE-KU (University of Copenhagen) expedition to the recently abandoned Camp Century site measured firn density to a depth of 15.0 m (Kovacs, 1970;Dansgaard et al., 1973).No temperature measurements appear to have been collected on this expedition.KU returned to resurvey the Camp Century borehole three more times, in 1977times, in , 1986times, in and 1989times, in . In 1977, firn density was measured to a depth of 100.0 m, and firn temperature was measured at 10 m depth (Clausen and Hammer, 1988).In 1986, firn density was measured to a depth of 12.0 m, and the 10 m firn temperature was again measured (Gundestrup et al., 1987).In 1989, the borehole temperature profile was measured from between 70 m and 270 m depth (Gundestrup et al., 1993).The upper depth at which this temperature profile was started (70 m in 1989) was taken as equivalent to the upper depth of the original temperature profile (50 m in 1966).Firn densities from these three KU resurveys (1977, 1986 and 1989) are not readily available at present, as the handwritten logbooks are currently in long-term storage (personal communication, Jørgen Peder Steffensen).
The NASA Program for Arctic Regional Climate Assessment renewed United States glaciological research at Camp Century.In 1995, the University of Colorado established an automatic weather station at "GITS", 5.2 km SE of Camp Century (Steffen et al., 1996).This automatic weather station logged hourly firn temperature measurements with thermocouples installed to a variety of depths up to 10.0 m (Sampson, 2009).Unfortunately, these instruments only operated from May 1996to January 1997, April 2001to June 2002, and March 2006to April 2007.These measurements are subject to significant noise (Sampson, 2009;Vandecrux et al., 2020a).We remove the most prominent peaks using variance filters as detailed in Vandecrux et al. (2020a).Residual noise (within ±1 °C) remains but measurements can nevertheless be used to evaluate the general performance of our firn model.Snow density has been periodically measured in annual snow pits during service visits to the GITS station (1998, 2006, 2008 and 2012).In 1996, the Ohio State University measured two firn density profiles, to depths of 21.8 and 120.5 m, at GITS (Mosely-Thompson et al., 2001).Firn temperature measurements do not appear to be readily available from the expedition.Annual vertical strain over the 0-21 m depth interval was measured by Ohio State University at GITS between June 1995 and May 1996 (Hamilton and Whillans, 2000).These vertical strain rate data, which was used to correct GPS-observed vertical velocity, unfortunately no longer appear readily available.
In 2008, the United States National Science Foundation (NSF), in collaboration with the USACE, initiated a series of overland traverses from Thule as a more cost-effective way of re-supplying Summit Station with fuel and other heavy cargo (Lever et al., 2016).The NSF overland traverse route transits some distance south of Camp Century.NSF overland traverses were undertaken in 2008, 2010, 2011, 2012, 2014and 2016 (personal communication, Jennifer Mercer).USACE researchers have participated in some of these traverses to undertake opportunistic glaciological observations.Annual snow density was measured at sites "B" and "I", approximately 6 km south of Camp Century, on the 2011 traverse (Hawley et al., 2014;Wong et al., 2015).No firn temperature measurements are readily available from these traverses.A temperature of −23 °C was measured at 9.5 m depth at Camp Century during the 2013 traverse (Polashenski et al., 2014).Firn density was measured to a depth of 10 m on the 2013 traverse, and to 2 m depth on the 2014 traverse (Polashenski et al., 2014), but these data do not appear readily available.In 2010, KU returned to Camp Century and measured firn density to a depth of 35.0 m (Buchardt et al., 2012).No firn temperature measurements were conducted during that expedition.
In 2017, the Geological Survey of Denmark and Greenland (GEUS) installed an AWS at Camp Century ("CEN"), with a similar design as the PROMICE AWS (Ahlstrøm et al., 2008) and that continuously measures a variety of climate variables, including firn temperatures to an installed depth of 10.0 m.This station is supplemented by a deeper thermistor station ("CEN-THM"), which continuously measures firn temperatures to installed depths of 62.0 and 73.0 m (Colgan et al., 2018).These deeper thermistors provide the first complete depth profile of firn temperatures within the uppermost 50 m of firn at Camp Century; above the upper limit of the original and re-surveyed borehole temperature profiles (Weertman, 1968;Gundestrup et al., 1993).GEUS also measured firn density to 73.0 m depth at Camp Century proper, within the debris field of the main tunnel network, as well as to a depth of 63.0 m at site 1.5 km upwind to the Southeast, where the firn density profile was undisturbed (Karlsson et al., 2019).
Finally, GEUS installed vertical strain gauges to measure daily firn compaction over the 0-5 m, 0-20 m and 1.4-62.3m depth ranges from August 2017 ("CEN-COM"; Colgan et al., 2018).The firn compaction stations are designed after Arthern et al. (2010).Each instrument consists of a line which is, on one end weighted and lowered into a vertical borehole and on the other end connected to a spring-loaded reel equipped with a potentiometer.As the snow and firn surrounding the borehole compact vertically, the borehole shortens, and the line is reeled gradually.This system allows to track the thickness of a firn portion through time.The daily compaction rates are calculated from borehole length change on daily resolution.The borehole remains "open"or air filledaround the potentiometer wire during these measurements.After approximately 2 years, however, deformational closure of the borehole, and/or exceeding the draw limit of the potentiometer, puts an end to the measurement.Due to the slow change rate of borehole length and the capacity of the datalogger to transmit borehole length at a ∼2 mm resolution, the compaction rates present spikes that are smoothed using a one-week moving window mean.
While we have endeavored to be comprehensive in this compilation of historical firn observations at Camp Century, we invite ongoing correspondence regarding any corrections or omissions, in order to improve the readily available in situ data for subsequent simulations of the firn at Camp Century (Table 1).

Climate Forcing (1966-2100)
To describe the surface climate at Camp Century from 1966 to 2100, we use simulations from the RACMO2.3p2regional climate model (Noël et al., 2019) and the CanESM2 earth system model (Chylek et al., 2011).The Regional Atmospheric Climate Model (RACMO2.3p2) is a high resolution (5.5 km) model forced by three different reanalysis dataset: ERA-40 over the period 1958-1978(Uppala et al., 2005)), ERA-Interim over 1979-2018 (Dee et al., 2011) and ERA5 during 2019 (Hersbach et al., 2020).The Canadian Earth System Model (CanESM2) combines the CanAM4 atmosphere models, CanOM4 ocean model, CanSIM1 sea-ice model and couples them to terrestrial and ocean carbon models (CTEM1 and CMOC1.2).Here, we use the model output produced for the CMIP5 experiment (Taylor et al., 2012) where atmospheric composition (including CO 2 concentration), solar forcing, aerosols and land use are prescribed.Over the 1966-2012 period, the historical CanESM2 output is used while for the 2013-2100 period, three CanESM2 simulations corresponding to RCP2.6, RCP4.5 and RCP8.5 are employed.CanESM2 produces five ensemble members (r1i1p1, r2i1p1, r3i1p1, r4i1p1 and r5i1p1) under each RCP scenario to explore the climatic variability within a given scenario.As our interest is in exploring variability between scenarios, we only use the r1i1p1 member of each scenario.We henceforth refer to both RACMO2.3p2and CanESM2 as "climate models".
The temporal resolution of climate model simulations can present challenges.Firn models best simulate melt and meltwater dynamics when using an hourly time step.Daily data from CanESM2 are consequently downscaled into hourly time steps using the following approach.For air temperature, relative humidity and wind speed, daily mean, maximum and minimum values are used to scale a daily sinusoidal function.This sinusoidal function oscillates around the daily mean given by CanESM2 and is prescribed to reach maximum value at 2:30 pm local time for air temperature, and at 4:30 pm local time for relative humidity and wind speed.Conversely, these functions reach their corresponding minima at 2:30 am and 4:30 am, respectively, local time.Illustration of these functions, as well as an assessment of their capacity to reproduce hourly values from daily values are available in Supplementary Figure S1.RACMO2.3p2data comes at a three-hourly temporal resolution and is linearly interpolated to hourly time steps.
Due to an imperfect representation of reality, climate model simulations differ from in situ observations.We therefore adjust climate model outputs to match in situ observations from the GITS and CEN stations.Both stations monitor on hourly intervals the air temperature, wind speed and humidity above the surface, as well as air pressure, downwelling and upwelling shortwave radiative fluxes.The GITS station monitors air temperature, humidity and wind speed at two heights above the surface.We preferably use the upper level, which is less prone to burial under accumulating snow, unless unavailable.The CEN station measures downwelling and upwelling longwave radiative fluxes unlike the GITS station.The GITS station provides data from 1995 to 2016 while CEN station covers 2017 to 2019.Air temperature, humidity, downward shortwave radiations, air pressure, wind speed and downward longwave radiation from the climate models are adjusted on a monthly basis so that the monthly mean and standard deviation of the adjusted values match the monthly means and standard deviations observed at GITS and CEN.A total of 8 years of observational data are available for air temperature, 9 years for shortwave radiation, 7 years for humidity, wind speed and air pressure, and about 2 years of measurements are available to adjust downward longwave radiation.We acknowledge that these are limited observations for the adjustment of our 54-135 yearlong climate model outputs.Nevertheless, the available weather observations span over a large range of values due to seasonal variations (Figure 2).These variations allow us to evaluate and adjust the climate model outputs for various climatic settings, which are, to a large extent, also representative of past and future climatic conditions at the site.The comparison of the adjusted simulations with observations are given in Figure 2. The comparison of the unadjusted climate model outputs to observations is presented in the Supplementary Figure S2.Different statistics apply for each RCP scenario because CanESM2 outputs differ after 2012.
For the upward shortwave radiative flux, we need an estimation of the surface albedo that can be extended into the future scenarios.Following previous work using local air temperature as a proxy for surface albedo (van As et al., 2013), we build an air-temperatureand solar-zenith-angle-dependent parameterization of surface albedo, α(T a , SZA), that we fit to Camp Century's nearest daily MODIS-observed albedo (Box et al., 2017).The parameterization reads as: α(T a , SZA) 422.8 − Ta 0.2 + Ta 2 53.2 -Ta 3 41384.4+ SZA 191.1 and is illustrated in Supplementary Figure S3.Eventually, we estimate upward shortwave radiation for each climate model by multiplying the adjusted downward shortwave radiation and the parameterized surface albedo.
Snowfall is a critical variable in firn simulations, as the vertical advection of mass within the model domain is ultimately controlled by long-term average snowfall rate.We also adjust the snowfall of climate model simulations to in situ, firn-corederived accumulation records: two cores covering 1966-1974 (Clausen and Hammer, 1988), two cores describing 1966-1998(Mosley-Thompson et al., 2001) and one core covering 1966-2006 (Buchardt et al., 2012).To minimize the squared difference from the median of these multiple annual accumulation observations, the hourly snowfall values are increased by 7.7% for RACMO2.3p2and by 12.8, 10.7 and 11.1% for the RCP2.6, 4.5 and 8.5 CanESM2 simulations, respectively.After the adjustment of air temperatures in the climate models, we reclassify the precipitation phase in each of the four datasets with a liquid-to-solid temperature threshold at 0 o C.

Surface Energy Budget, Mass Balance and Firn Model
For each adjusted climate model dataset, the surface energy budget is calculated for every hourly time step as the sum of the downward and upward shortwave radiation, downward and At 10 m Gundestrup et al. (1987Gundestrup et al. ( ) 1966 At 50 m Weertman (1968), Gundestrup et al. (1993Gundestrup et al. ( ) 1996Gundestrup et al. ( -1997 10 m long thermocouple string at GITS, GC-Net Steffen et al. (1996Steffen et al. ( ) 2001Steffen et al. ( -2002 Depth  1988;unstable: Paulson, 1970;Dyer, 1974) and using the surface roughness length parametrizations of Andreas (1987) and Lefebre et al. (2003).The conductive heat flux is calculated by our firn model.The GEUS firn model is a multilayer snow and firn model (Langen et al., 2017;Vandecrux et al., 2018;2020a) where each of the 100 layers is composed of snow, ice and water compartments.During snowfall, new snow is added to the top of the model column with a density of 315 kg m −3 after Fausto et al. ( 2018).When melting is prescribed at the surface, snow and ice are shifted to the liquid water compartment of the first layer.The meltwater is then able to move vertically according to the parametrization of Darcy's law developed by Hirashima et al. (2010).If liquid water percolates into a subfreezing layer, the cold content of the layer is used to refreeze the percolating water and the mass of refrozen water is added to the ice compartment of that layer.The model accounts for thermal diffusion using the thermal conductivity and the heat capacity of snow from Yen (1981).Each hourly time step, the firn density is updated for compaction according to the weight of overlying snow and firn (Vionnet et al., 2012).The grain size is also calculated using the formulation from Brun (1989).The snow and firn model is initiated with the observed density profile from 1966.In the absence of an observed temperature profile from 1966, we use the temperature profile from Site 2, surveyed in 1957, located 130 km away from Camp Century but presenting identical average accumulation and temperature as Camp Century (Spencer et al., 2001).The model is used to track the burial of the upper limit of the debris field as detected at 32 m depth in 2017 by Karlsson et al. (2019) as well as to track the maximal meltwater percolation depth each year.

Climatology, Surface Energy and Mass Budget
Following the adjustment of climate model outputs, we obtain an overview of both past and future climate at Camp Century over the 1966-2100 period.This adjusted climate forcing, which is used to force the GEUS firn model, are displayed in Figure 3 and their mean, standard deviation and decadal trends over the 1966-2019 and 2020-2100 periods are reported in the Supplementary Table S1.Clear positive trends in air temperature, 0.38-0.54°C decade −1 , are simulated by all models for the 1966-2019 period.This warming is predicted to accelerate in the RCP8.5 scenario, to reach a plateau in RCP4.5 and to slow down in RCP2.6, with trends of 0.76, 0.33 and 0.17 °C decade −1 over the 2020-2100 period.These trends correspond to increases in mean annual air temperature of 6.0, 2.6 and 1.4 °C over the next 80 years.Other positive trends in relative humidity and downward longwave radiative fluxes indicate a warmer, more humid atmosphere at Camp Century in the future.These trends are more pronounced in the RCP8.5 scenario but rarely reach statistical significance (p-value < 0.01) due to high interannual variation (Supplementary Table S1).Wind speed exhibits no statistically significant trends over the whole period.Downward shortwave radiative fluxes are relatively stable over the historical period and decrease slightly during 2020-2100.The temperaturedependent parametrization of surface albedo predicts a perceptible decrease in surface albedo in the CanESM2 RCP8.5 simulation after 2050.
We use the adjusted climate variables to calculate the surface energy budget at Camp Century during the 1966-2100 period (Figure 4).For each energy flux, the mean, standard deviation and decadal trend over the 1966-2019 and 2020-2100 periods are reported in the Supplementary Table S1.The sensible heat flux decreases slightly but remains positive in the four model runs over the 1966-2019 period.This decrease of sensible heat flux indicates that in spite of warming temperatures in all scenarios (Figure 3), the temperature gradient between the air and the surface decreased over the entire study period.Indeed, the surface temperature rose by ∼2.5 °C between 1966 and 2019 and by another 2-6 °C, depending on the scenario, between 2020 and 2100 (Figure 4).The annual average latent heat flux is negative and decreasing for all model outputs indicating an increasing sublimation, also due to increasing surface temperatures (Figure 4).The mean conductive heat flux from the surface to the deeper firn remains stable with lowest mean value in the RCP8.5 scenario.This is due to the simultaneous warming of the surface and underlying firn, which combined, lower the temperature gradients in the near surface firn and reduce heat transfer from the surface to the firn.Rain is absent from the adjusted model runs in the 20 th century and appears only in the RCP4.5 and RCP8.5 scenarios after ∼2065.The magnitude of these events and the energy that they bring to the surface remain marginal (Figure 4).The net shortwave radiative flux is the dominant source of energy for the surface.It increases most in the RCP8.5 scenario due to a stronger decline of albedo in that model output.The net longwave radiation is the main energy sink for the snow surface.Its magnitude decreases through the entire study period, indicating that the increasing emitted longwave radiation does not compensate for the increasing incoming longwave radiation.The sum of these fluxes results either in a change of the snow temperature or, in the summer, to snow melt.Melt increases at a rate of 1.7-3.8mm water equivalent (w.e.) per decade during the historical period.Over the 2020-2100 period, this increase either continues at similar rates in RCP2.6 and RCP4.5 scenarios or even steepens to 10.6 mm w. e. decade −1 in RCP8.5.We note that the adjusted RACMO2.3p2output captures melt peaks in 2010 and 2012 while the adjusted CanESM2 data do not.
The firn model, in combination with any of the four forcing datasets, do not calculate meltwater runoff in either the historical or future period.Consequently, sublimation is the only negative component of the surface mass budget.Following the increase of surface temperature, sublimation is expected to remove increasing amounts of mass from the surface (Figure 5B) and will represent ∼10% of snowfall by the end of the century.Snowfall is relatively stable in all the adjusted model outputs apart from RCP8.5 which shows a slight decrease (Figure 5A).These mass fluxes result in a surface mass balance that does not show any statistically significant trend apart for the RCP8.5 scenario where surface mass balance decreases by 7 mm w. e. decade −1 over the 1966-2100 period.
The last estimation of surface mass balance at Camp Century dates back to Colgan et al. (2016) (Figure 6).Colgan et al. (2016) used MAR regional climate model forced by the same CanESM2 outputs as here, but without adjusting any climate variable other than snowfall for potential biases.Snowfall rates employed by both studies are therefore very similar (Figure 6).In contrast, melt rates presented by Colgan et al. (2016) during the historical period are up to an order of magnitude greater than the melt rates calculated in our study.We attribute this melt-rate discrepancy to the positive mean bias of unadjusted CanESM2 air temperatures at Camp Century: between +3.7 °C and +5.1 °C depending on the RCP scenario and evaluation AWS (see Supplementary Figure S2).While our present study corrects this air temperatures bias based on observations, this correction was not performed for the MAR-derived melt rates employed by Colgan et al. (2016).Unadjusted CanESM2 air temperatures depict a melt season at Camp Century that is longer, and warmer, than is realistic.Differences in melt rates also explain the sharp contrast in runoff components between the studies: none here, while Colgan et al. (2016) predicted that runoff would exceed annual snowfall by the end of the century (Figure 6).We also note that different runoff parameterizations are used in the two studies, resulting in differing fractional runoffs for a given melt rate.

Evolution and Evaluation of Firn Temperature
The capacity of the model to simulate the meltwater percolation and the firn density greatly depends on how accurately it can simulate  firn temperature (Figure 7).While all models reproduce the warming trend visible in 10 m firn temperature observations, they underestimate punctual measurements by up to 0.5 °C.Following the increased air temperature and downward longwave radiation present in the surface forcing dataset, the firn undergoes warming during the 21 st century at a rate of 0.2-0.3°C decade −1 over 1966-2019.During 2020-2100, the 10 m firn temperature increases at a rate of 0.2, 0.3, 0.5 °C decade −1 and reaches −21.3, −20.0 and −18.6 °C by 2100 in the RCP2.6, 4.5 and 8.5 scenarios respectively.The comparison on hourly time steps against GITS and CEN_THM firn temperature observations confirms the cold bias of all models with an average bias of −0.2 °C at GITS and up to −0.9 °C at CEN_THM.Yet the RMSE is within one degree at CEN_THM, with a R 2 above 0.95 while it reaches 2.5 °C at GITS, with R 2 of about 0.4 for all models, probably due to low signal-to-noize ratio within the GITS dataset.The cold bias is nevertheless less pronounced during summer season.

Evolution and Evaluation of Firn Density
In spite of a firn warming at a rate of 0.2-0.3°C decade −1 during the historical period confirmed by punctual firn temperature observations (Figure 7), the firn density has remained relatively stable.Indeed, the observed firn densities in 1966, 1996, 2010 and 2017 all fall within 40 kg m −3 of each other (Supplementary Figure S5).This stability of firn densities indicate that the densification rates are insensitive to changes in firn temperatures for these specific firn conditions.This is captured accurately by the firn model since these punctual observations match simulated firn density profiles down to 30 m depth (Figure 8).Nevertheless, although our simulations match relatively well the observed firn density in 1996, they all underestimate the firn densities below 30 m in 2011 and 2017 (Figures 8F-H).We attribute this discrepancy to the use of the compaction parametrization from Vionnet et al. (2012), that was designed for seasonal snow and may underestimate compaction rates below 30 m.We also note that single firn density profile observations are affected by the spatial heterogeneity of the firn (e.g., Steffensen et al., 1996), which was not evaluated at Camp Century.

Evaluation of Firn Compaction Rates
Another key process that determines the burial rate of the military waste at Camp Century is the gradual compaction of the firn.A model that overestimates firn compaction rates will  (1976,2011) and by thermistor strings (2002,2006,2017).(F-H) Punctual comparison of the simulated (colored line) and observed (black diamonds) temperature profile.
underestimate the depth of the debris field.The simulated compaction rates depend on simulated surface accumulation (Figure 5C), firn temperatures and densities (Figures 7, 8) but also on the choice of the compaction formulation (Vionnet et al., 2012).They are evaluated using the CEN-COM instrument, for which the top and bottom depth of the boreholes are also tracked within the firn model (Figure 9).The deepest borehole, spanning over 1.4-62.3m depth, shortens slower in the simulations compared to the observations (Figure 9A).Consequently the simulated compaction rates are 24-38% lower in the simulations compared to the observation (Figure 9E).When compared to the instruments installed over 0-20 m and 0-5 m depth (Figures 9B,C), models show more differences in performance.The models have mean biases ranging from 2 to 25% at the 0-20 m instruments (Figures 9B,F) and from −30 to 7% for the 0-5 m instrument (Figures 9C,G).Lastly, the models have a mean bias within ±9% when evaluated at the shortest boreholes, installed at 0-4.9 m depth (Figures 9D,H).The larger mean bias at the longest instrument and smaller mean bias at the shortest instrument support the hypothesis that the compaction scheme from Vionnet et al. (2012) is more accurate at describing the nearsurface snow and less accurate for deeper firn.Nevertheless, the R 2 between simulated and observed compaction rates ranges between 0.65 and 0.93 depending on the instrument-model combination.This confirms that the model captures the seasonal variations of the firn densification for all the instruments.The spread between models and the high temporal variability in the three shorter instruments also reveal the effect of snowfall, which magnitude and timing are different in each simulation.
Beyond the simple evaluation of our model's capacity to simulate firn compaction and burial rates, our in situ compaction observations highlight that the relevant depth range for in situ monitoring of firn compaction is significantly deeper than typically assumed.Previous studies seeking to measure firn compactionin order to constrain local mass balancehave only measured firn compaction within the top c. 20 m of the firn column (Hamilton and Whillans, 2000;Burgess and Sharp, 2008).Our in situ measurements show that, at least in 2018, the average compaction rate was 25% higher in the 62.3 m deep borehole (−1.31 mm day −1 ) than in the 20 m deep borehole (−1.05 mm day −1 ).This highlights that significant compaction occurs below the typical 20 m monitoring depth.While this is consistent with theory, our firn model simulates the seasonality of compaction rates in the 62.3 m borehole but underestimates their magnitude.This further underlines the value of the deep firn compaction observations for calibrating firn compaction models and assessing local mass balance via in situ or remote sensing approaches.To our knowledge, the CEN_COM instrument is the deepest installation of a compaction instrument to date on the Greenland ice sheet.

Depth of the 1966 and 2012 Horizons and Meltwater Percolation
Considering the satisfactory performance of the firn model in terms of firn temperature, density and compaction rates we track the 1966 and 2012 summer surfaces in the model outputs as they get buried under new snow accumulating at the surface and compacting with time.The depth of these horizons can be checked in 2017 when a radar campaign located the 1966 summer surface at 32 ± 2 m depth and the 2012 summer surface at 7 ± 0.5 m depth.Simulations tend to underestimate the depth of these horizons: respectively 28 ± 1 m and 4 ± 0.5 m depth.Nevertheless we consider these results encouraging and indicative of the suitability of the surface forcing and firn model set up for the tracking of the horizons in the future.
When tracking the top of the debris layer, observed at 32 m in 2017, over the next decades, the three RCP scenarios are consistent until the 2060s when they start departing from each other (Figure 10A).In comparison to the low and medium emission scenarios (RCP 2.6 and 4.5) the RCP8.5 model run keeps the debris layer closest to the surface and its burial rate slows down but does not revert by the end of the simulation in 2100.The firn model suggests that the debris field upper horizon will continue to be buried by persistent net accumulation over the next eighty years under all RCP scenarios.The RCP2.6 scenario suggests that the debris field upper horizon will be located at 64 m depth in 2100, while the RCP8.5 scenario suggests an analogous depth of 58 m.This persistent burial reflects projections of positive surface mass balancedespite increasing meltunder all RCP scenarios to 2100.This estimate remains of the same order as the one from Colgan et al. (2016) who estimated main trench depth as 67 m in 2090 using unadjusted climate model outputs under RCP8.5 scenario.The maximum simulated depth of meltwater percolation by 2100 is about 1.1 m in RCP8.5 and 0.6 in RCP2.6 and RCP4.5 (Figure 10B).This is almost two orders of magnitude less than the depth of the debris field upper horizon projected in any scenario at that time.We therefore find extremely unlikely that meltwater would interact with the subsurface debris field at Camp Century before 2100 under all RCP scenarios.This conclusion is dependent on assumptions embedded in the parameterization of physical processes within the GEUS firn model, as well as assumptions in prescribed climate forcing.For instance, the GEUS firn model parameterizes the vertical percolation of meltwater with a Darcy flow formulation.The firn meltwater Retention Model InterComparison Project (RetMIP) suggests that other parameterizations (e.g., bucket schemes or models including heterogenous percolation) can result into meltwater percolation as much as five time deeper than the GEUS firn model under the same climate forcing (Vandecrux et al., 2020b).Yet a five-fold increase in meltwater infiltration depth would still be insufficient to bring meltwater at the depth of the debris field.The impact of the choice of forcing is visible when comparing our results with the ones from Colgan et al. (2016).Indeed their unadjusted CanESM2 forcing predicted net ablation, and therefore upward movement of the debris field by the end of the century.Additionally, we note that although the adjusted CanESM2 and RACMO2.3p2forcing datasets presented similar melt rates over the historical period, the percolation depth in the RACMO2.3p2-forcedmodel output reaches 0.9 m in 2012, twice as much as the maximum depth of simulated by CanESM2driven models that same year.We hypothesize that the RACMO2.3p2output (available at three-hourly values) is able to capture more intense melt rates, and therefore meltwater percolation events, than the daily CanESM2 output, even after our temporal downscaling.Nevertheless, we consider this difference between models insufficient to influence the possibility of meltwater interaction with the debris field before the end of the century.

SUMMARY REMARKS
This study explores past and future evolution of the firn at Camp Century in order to assess the possibility of meltwater interacting with the subsurface debris field.We employ the GEUS firn model to simulate firn density and temperature within the uppermost 80 m of firn over the 1966 to 2100 period.The model is forced by regional and global climate model simulations that were adjusted to automatic weather station data and to net accumulation reconstructed from five firn cores.
Under this prescribed climate forcing, the GEUS firn model reproduces, within acceptable uncertainty and biases, a variety of in situ observations: firn density measured in three firn cores, firn temperature measured at depths ranging from 0 to 73 m, firn compaction rate measured over multiple depth ranges, and depth of the 1966 horizon as observed in 2017.This robust evaluation of firn model performance during the historical period  lends confidence to the firn model output during the future period (2020-2100).
The firn model suggests that firn temperatures will continue to warm at Camp Century under the full range of IPCC Representative Concentration Pathways (RCPs).Under the RCP2.6 scenario (respectively RCP8.5), the firn at 10 m depth will warm from −24.0 °C in 1966 to −21.3 °C (resp.−18.6 °C) in 2100.Under these scenarios, Camp Century will transition from a dry snow zone, where negligible melting occurs, to a percolation zone, where melt occurs and refreezes within the annual snowpack.This will be accompanied by the emergence of discrete annual ice layers within the firn.As a result of persistent positive mass balance, the debris field is predicted to be below 58 m of firn and still being further buried at the end of the century in all RCP scenarios.In combination with the absence of meltwater percolation below 1.1 m in all simulations, it is therefore extremely unlikely that meltwater interacts with military waste within this century.
While the groundbreaking glaciology and climatology research undertaken at Camp Century in the 1960s provides an almost unparalleled scientific heritage amongst Greenland ice sheet research sites, the military history of Camp Century now gives the site unanticipated social significance in light of climate change.Our study represents the first assessment from the Camp Century Climate Monitoring Program on the possibility of meltwater interacting with the subsurface debris field at Camp Century over the next eighty years.

FIGURE 1 |
FIGURE 1 | Location of Camp Century.Ice sheet is represented by the gray area, 1000 m elevation contours by thick white lines and 250 m contours by thin white lines.
correction as in Vandecrux et al. (2020b) 2006-2007 2017 (ongoing) 73 m long thermistor string, Colgan et al. (2018) Firn compaction 2017 (ongoing) Three instruments operating from August 2017 until May 2019 at 1.4-62.3m, 0-20 m, 0-5 m depth, one instrument operating since May 2019 at 0-4.9 m depth, Colgan et al. (2018) upward longwave radiation, latent and sensible turbulent heat fluxes, and conductive heat flux through the top layer of snow.For each time step, we find iteratively the surface temperature that satisfies the closure of the surface energy balance.When the sum of energy fluxes is positive and the surface snow is at melting point, the surface temperature remains at 0 °C and the excess energy from the surface energy budget is allocated to the melting of snow and ice.Downward shortwave and longwave radiative fluxes are directly taken from the adjusted climate model data.Upward shortwave radiation is calculated using our parametrization for albedo (Section 2.2) and upward longwave radiation is calculated from the surface temperature using the Stefan-Boltzmann law and a surface emissivity of 0.98.The latent and sensible heat fluxes are calculated from air and surface temperatures, humidity and wind speed after van As et al. (2005) as implemented in Vandecrux et al. (2018), Vandecrux et al. (2020a).The sensible and latent heat fluxes calculation is based on the Monin-Obukhov similarity theory and accounts for the boundary layer stratification (stable: Holtslag and Bruin,

FIGURE 2 |
FIGURE 2 | Comparison of hourly meteorological fields from adjusted climate models with automatic weather station observations.

FIGURE 4 |
FIGURE 4 | Yearly averages of the heat fluxes contributing to the surface energy budget (positive when warming the surface) at Camp Century over the 1966-2100 period: Sensible heat flux (A), latent heat flux (B), conductive heat flux from subsurface (C), rain heat flux (D), net shortwave (E) and longwave (F) radiative fluxes and the resulting surface temperature (G), and total annual meltwater production (H).

FIGURE |
FIGURE |Yearly total of the surface mass balance components at Camp Century over the 1966-2100 period: snowfall (A) and net sublimation (B).Note that negative values for sublimation indicate mass removal from the surface.No runoff was calculated.Surface mass balance (C) in the adjusted models and observed of five firn cores.

FIGURE 6 |
FIGURE 6 | Snowfall (A), surface melt (B), runoff (C) and surface mass balance (D) used in Colgan et al. (2016): calculated by MAR when forced by unadjusted CanESM2 (blue lines), or NorESM1 (green line) outputs; and in this study: calculated by the GEUS model forced by adjusted CanESM2.

FIGURE 8 |
FIGURE 8 | (A-D) Simulated density in the top 40 m of firn.(E-I) Punctual comparison of simulated (colored line) and observed (black line) firn density.

FIGURE 9 |
FIGURE 9 | (A-D) Observed (black line) and simulated (colored lines) borehole length evolution.(E-H) Observed (black line) and simulated (colored lines) daily firn compaction rate in each borehole.

TABLE 1 |
Overview of the firn measurements at Camp Century used for the evaluation of our firn model.