New large subglacial lake in Princess Elizabeth Land, East Antarctica, detected by airborne geophysical observations

Knowledge of subglacial lakes is important for understanding the stability of the Antarctica Ice Sheet (AIS) and its contribution to the global sea-level change. We designed an intensified airborne campaign to collect geophysical data in Princess Elizabeth Land (PEL), East Antarctica, during the 2015-2019 CHINARE expeditions. We developed an innovative method to build a set of evidence of a newly detected subglacial lake, Lake Zhongshan. 15 Adaptive RES data analysis allowed us to detect the lake surface and extent. We quantified the lake depth and volume via gravity modeling. Another dataset collected at Lake Vostok provided the ground truth. The results revealed that Lake Zhongshan, located at 73°26'53"S, 80°30'39"E and ~3,603 m below surface, has an area of 3281 km, making it the only one in PEL and the fifth largest in Antarctica. These findings are important for understanding subglacial hydrodynamics in PEL, as well as the stability of the AIS. 20


Introduction
Antarctic subglacial lakes are reservoirs of water gathered beneath the Antarctic Ice Sheet (AIS) (Christianson et al., 2012;Napoleoni et al., 2020). Large subglacial lakes can affect the overlaying internal ice sheet structure, rheology of the ice body (Siegert. et al., 2000;Pattyn et al., 2018), and modulate the velocities of ice streams and outlet glaciers (Robin, 2008). Some subglacial lakes could be connected over large distances through a hierarchical hydrological 25 chain (Siegert et al., 2016;Wingham et al., 2006). The growing evidence for a well-connected subglacial drainage system in Antarctica indicates that the subglacial hydrologic system is indispensable to understanding the mass balance of AIS and its contribution to global sea-level change (Jordan et al., 2010;Robin, 2008). Although the methodology of AIS mass balance has been widely studied Rignot et al., 2019;IMBIE, 2018), the role of the subglacial hydrological system in this ice sheet mass change process is poorly understood (IPCC, 2019). 30 Therefore, accurate detection of subglacial lakes and investigation of their hydrological connectivity are essential to learning spatial and temporal evolution of the overlaying ice sheet dynamics.
Subglacial lakes are usually recognized as a flat interface with high specularity in radargrams from radio echo sounding (RES) data (Oswald and Robin, 1973;Wright and Siegert, 2012). Until the early 1990s, RES data were the https://doi.org/10.5194/tc-2021-332 Preprint. Discussion started: 3 December 2021 c Author(s) 2021. CC BY 4.0 License. only source used to detect subglacial lakes (Siegert, 2018). Satellite altimetry data, including ERS-1, ICESat, CryoSat-35 2, and ICESat-2, provide another important source for detecting "active" subglacial lakes since the elevation changes caused by subglacial hydrologic activities can be detected by the time-series data from satellite altimetry (Malczyk et al., 2020;Siegfried and Fricker, 2018). However, this method cannot be used to detect "definite" lakes (Rivera et al., 2015). Therefore, RES observations are still the most reliable and irreplaceable data for the detection of both active and definite subglacial lakes. Additionally, airborne gravity and seismic sounding observations have been applied to 40 determine lake bathymetries (Studinger et al., 2004).
With great efforts in last five decades, the inventory of Antarctic subglacial lakes has been updated in a fourth edition (Wright and Siegert, 2012), and the tally of known discrete subglacial lakes has reached more than 400 (Siegert et al., 2016). However, for a long time there has been a lack of geophysical observations in the vast sector of Princess Elizabeth Land (PEL), East Antarctica, which lies between the Amery Ice Shelf and Wilhelm II Land (Cui et al., 45 2020b). Whether there are subglacial lakes in PEL is unknown (Fretwell et al., 2013). Under the hypothesis that ice surface expressions are closely related to subglacial features, Jamieson et al. (2016) inferred that a large subglacial lake may exist in the interior of PEL (Fig. 1b) as an elongate, extensive, and relatively featureless zone of the ice sheet surface appeared in MOA and RADARSAT imagery. This statement has not yet been verified by independent and direct geophysical observations. 50 We designed and carried out a geophysical data campaign during the Chinese National Antarctic Research Expedition (CHINARE) from 2015 to 2019, which aimed at filling the geophysical data gap in PEL under the auspices of the ICECAP2 program (ICECAP2, 2021), but targeted the potential subglacial lake region with intensified observations (Fig. 1a). A fixed-wing aircraft, Snow Eagle 601, collected a suite of repeated geophysical data, including RES, gravity, and magnetic data in PEL (Cui et al., 2020a). In this paper, RES and gravity data of 21 survey lines in 55 an extended region of the potential lake are used for subglacial lake detection. The RES data in an area of seven survey lines (black line segments in Fig. 1b) show a strong water signature of basal reflectivity and prove the existence of this subglacial lake, now named Lake Zhongshan. Furthermore, we developed methods for the in-depth analysis of RES data and gravity inversion modeling which are adaptive to our data acquisition system and the local subglacial environment in PEL. We also collected another dataset from Lake Vostok using the same sensor system to provide 60 the ground truth from a known Antarctic subglacial lake to validate our results for Lake Zhongshan. The results of this paper reveal the geometric attributes of the new lake, including location, elevation, boundary, area, volume, and bathymetry, among others. In addition, the uncertainty of the estimated and modeled parameters is analyzed. These  The average reflectivity difference between segment AA' and surrounding areas is used to determine the high confidence segment for the lake area.

Data
Snow Eagle 601 has successfully collected airborne RES data along flight lines of more than 175,000 km during the PEL geophysical data campaign from 2015 to 2019 (Cui et al., 2020a). The penetration capability of these RES 80 instruments has been confirmed by the detection of Lake Vostok with an ice thickness of ~4,200 m (Cui et al., 2018). lines. The average flight altitude was ~600 m above the ice surface to guarantee sufficient penetration capability and avoid signal saturation for a high gain radar channel (Cui et al., 2020a). The radar dataset was sampled every ~17 m along-track, and the depth resolution is ~5.6 m in ice. Raw data collected from the survey flights were converted into 90 readable datasets and georeferenced by using onboard near real-time GNSS positions through the Environment for Linked Streams Acquisition (ELSA) system developed by the University of Texas Institute for Geophysics (UTIG) (Ng et al., 2020). The data quality was carefully controlled through both high and low gain channels.
The gravity data were acquired on the same survey flights as the RES data by a GT-2A airborne gravimeter (Cui et al., 2018). This three-axis stable scalar gravimeter can perform high-precision gravity data collection under turbulent 95 air currents and varying flying heights. Four JAVAD Delta dual-frequency and carrier phase GNSS receivers were installed at four different positions on the aircraft (two wings, front and center of the fuselage), providing a positional accuracy of 25 cm under normal circumstances (Smoller et al., 2015). In addition, the connection between the GNSS antennas and the GT-2A system was boresight calibrated. Thus, the gravity data (and further the RES data) were georeferenced by the GNSS system in both time and space. The gravity data sampling frequency is 2 Hz. Free-air 100 gravity anomalies were filtered with a 200 s full wavelength filter, resulting in a ~8.3 km half-wavelength resolution for the flight speed of ~300 km/h (~83 m/s). A Kalman filter was applied in real-time to eliminate the effect of additional horizontal acceleration during flight.

Detection of subglacial lake surface using basal radar reflectivity 105
Detection of subglacial lakes is mainly realized by analyzing the radar reflectivity from the bed of the ice caps in RES radargram profiles (Fig. 1c) and further identifying the ice-water interface as the lake surface. To distinguish potential lakes from their surroundings, we rearrange the radar equation (Bogorodskiy et al., 1985;Peters et al., 2007) for estimating reflection R as follows: https://doi.org/10.5194/tc-2021-332 Preprint. Discussion started: 3 December 2021 c Author(s) 2021. CC BY 4.0 License.
where and denote the transmitted and received power, respectively, is the antenna gain, is the wavelength of the carrier frequency, is the transmission coefficient at the air-ice interface, is the reflection loss of the icebedrock interface, is the system gain, is the ice thickness, ℎ is the height of the aircraft above ground, and is the ice attenuation loss. The received power is obtained by integrating power across ±800 ns of the picked bed position (Humbert et al., 2018). The ice attenuation loss L A is integrated over the travel path between the surface h s 115 and the base h b (Matsuoka et al., 2010): where the A(z) denotes the local attenuation rate per unit path length in the vertical (z) direction. We quantify the A(z) following previous studies (e.g., Humbert et al., 2018;MacGregor et al., 2007;Matsuoka et al., 2010Matsuoka et al., , 2012Zirizzotti et al., 2010): 120 where 0 = 8.8541878176 × 10 −12 F/ m −1 and 0 = 1.25663706 × 10 −6 N/ A 2 are the free space dielectric permittivity and magnetic permeability, respectively. The conductivity ( ) is calculated by the Arrhenius model as follows: where E is the activation energy, k is the Boltzmann constant, and 0 is the reference conductivity for ice, which can be affected by the impurities. In this study, we set 0 at = 251 K (Matsuoka et al., 2012). We finally visually evaluate variations of the reflectivity that represent the returned power from the bottom of the ice-water and ice-rock interfaces. Any region with an abrupt increase in bed return power that exceeds a threshold ranging from 2.0 dB  to 7.7 dB (Zirizzotti et al., 2010) may be a safe indicator of a subglacial water body. 130

Lake depth estimation by gravity inversion
We use the GT-2A gravity data in the same map extent as the RES dataset to estimate the lake depth through gravity inversion. The Parker inversion method (Parker, 1972) implemented in the Geosoft GM-SYS 3D system is used to iteratively approach the lake bottom by minimizing the misfit between the calculated and observed gravity values.
In principle, the variation in a gravity field can be observed and presented as gravity anomaly ∆g that can be further 135 modeled from the density contrast across the interface between two layers ρ, depth to the interface h(x), among other parameters. Its Fourier transform is expressed as where G is the universal gravitational constant, k is the wave number, h(x) is positive downward, and z0 is the mean depth of the horizontal interface. The model domain is represented as three horizontal layers: a) a solid ice layer with 140 a density of 0.917 g/cm 3 , b) a lake water layer with a density of 1.028 g/cm 3 , and c) a rock layer with a density of 2.67 g/cm 3 . In a prior study, we showed that this three-layer model provided superior performance compared to a four-https://doi.org/10.5194/tc-2021-332 Preprint. Discussion started: 3 December 2021 c Author(s) 2021. CC BY 4.0 License. layer model, which includes an additional fresh sediment layer or treats the rock density as a variable . A forward gravity model is calculated using BedMachine Antarctica data (Morlighem et al., 2020) as the initial bed solution (Fig. 3c). We then calculate the direct current (DC) shift or difference between the modeled and observed 145 gravity at locations where the ice thickness is known from the RES data (Fig. 1b). At each location along the survey line, the DC shift exhibits the vertical spatial variation that reflects the varying water layer used to approach an optimal lake depth. We then produce a regular grid of the DC shift using a minimum curvature interpolation algorithm (Smith and Wesse, 1990, Fig. 4b). We then correct the observed gravity with the interpolated DC shift and use it as input to the inversion model to account for the original bias. Since we fill observational gaps between the survey lines with 150 the regionally corrected model data (Fig. 4c), the impact of the data gaps on the inversion process, especially at the edges of the survey region, is minimized. We allow an 8000 m wide horizontal transition to enable a smooth transition from the model to the observations so that all observations are preserved. The grid spacing is chosen as 1000 m. The lake depths are finally estimated by the gravity inversion method, which minimizes the misfit between the observed and modeled gravity. We also generate an error map for the lake bathymetry product (Fig. 4a). 155

Lake Zhongshanthe newly detected subglacial lake in PEL
Here we present the geophysical evidence of the newly detected subglacial lake, now named Lake Zhongshan, in PEL, East Antarctica.

Localization and attributes of Lake Zhongshan derived from RES observations
Initially, an extensive search for basal reflectivity anomalies based on Equations (1-4) was performed using all RES 160 data collected along the survey lines (Fig. 1a), especially along those around the optical lake boundary (Fig. 1b).
Finding of the vertical zone of the bedrock layer and determination of the ice-water interfaces were mainly carried out by visual interpretation of the features presented in the RES radargram profiles (Fig. 1c). Any relatively long and bright linear features, compared to that of the surrounding bedrock, are considered candidates for the ice-water interface. We then confirm a high confidence segment of a subglacial lake if its maximum reflectivity increase exceeds 165 7 dB and the average reflectivity increase is over 3 dB. For example, survey line segment AA' in Fig. 1b showed the longest bright stretch of 80 km in this region (Fig. 1c), with an average reflectivity increase of 4.8 dB and the maximum increase over 9.5 dB in the middle section (Fig. 1d). This is thus considered a high confidence segment (black line in Fig. 1b). Subsequently, six other high confidence segments were detected in the same way (Figs. 2 and 1b). The overall average reflectivity increase for all seven segments is 3.61.0 dB and their maximum reflectivity increases range from 170 7.2 to 10.2 dB, indicating a high level of confidence that the RES data along these segments provide the areal signature of a subglacial lake from a radargram point of view.  On the other hand, the combined RES survey lines of all three austral seasons in the region have an average spacing of 14 km between the lines. With three intensified survey lines, the area of the high confidence segments has 180 a densified average spacing of 7.6 km. Furthermore, we have a high radar sampling spacing of ~17 m along-track. Therefore, the seven high confidence segments are combined to form a "continuous" lake area enclosed by a boundary (red line in Fig. 1b) formed by linking the end points of the high confidence segments. Consequently, the lake's dimensions are determined to be 45,20028 m in length and 10,90028 m in width. The area within the lake boundary is 3281.5 km 2 . Thus, Lake Zhongshan is the fifth largest in Antarctica. We determined the center of the lake boundary 185 as the geographic location of the lake (73°26'53"S, 80°30'39"E) with an uncertainty of under 1 m. In addition, the bed elevation data show a basin topography in the detected lake area (Fig. 1b). Based on the ice cap surface model of REMA (Howat et al., 2019) and the RES data, we estimated that the lake surface lies ~700 m below sea level and the ice above the lake is ~3603 m thick. https://doi.org/10.5194/tc-2021-332 Preprint. Discussion started: 3 December 2021 c Author(s) 2021. CC BY 4.0 License.

Depth of Lake Zhongshan modeled from gravity data 190
The observed gravity anomalies along the survey lines within the area vary from -35 mGal in the RES-derived lake area to 55 mGal in the surroundings (Fig. 3a). The modeling results from the gravity data using Equation (5) present a subglacial lake whose boundary shape generally matches that of the RES lake boundary (Fig. 3b). The lake is asymmetric in the north-south direction and contains two non-flat lake basins with a water depth of 200-250 m. The northern basin is 50 m deeper and approximately double the area of the southern basin. The ridge that connects the 195 two basins is as shallow as ~150 m. We estimated a lake volume of 54 km 3 . Furthermore, the bed topography presented by the BedMachine Antarctica (Morlighem et al., 2020) was used as an initial interface between the ice and bedrock in the larger surrounding region of the lake for the gravity modeling process (Fig. 3c). Within the lake area, the gravityinferred bed (Fig. 3d) is on average ~200 m deeper than the initial bed The surrounding area outside of the optical lake boundary is not changed during the inversion process. This indicates that the collected airborne gravity data and 200 the gravity model fit well for the overall region, while the observed gravity anomalies within the lake area provide evidence of the existence of a water body.

210
We use the gravity misfit estimated as the difference between the modeled gravity and the observed gravity to quantify the uncertainty of the inversion (Fig. 4a). To translate the gravity misfit into an error in bed elevation, we performed a numerical experiment in which we calculated the differences between the gravity field in the region and

225
Along the survey lines of X61a_2016, X61a_2018, Y75a and Y81a (Fig. 1b), there are connecting areas among the lake segments where the reflectivity values are as high as on the ice-water interface (Figs. 2a-2d). However, the corresponding brightness features in the radargrams are of high roughness and are not as flat. They appear to represent shallow lake boundary areas of a wet sediment layer. In fact these areas are coincident with the extended shallow areas of the gravity-inferred lake outside of the RES lake boundary (red line) in Figure 2b. We used a 3D gravity modeling 230 method that does not include the sediment layer because of its resulting uncertainty. In future research, we will expand our gravity modeling with more components to determine the sediment layer of the lake.
We recommend using the RES lake boundary (red line in Fig. 1b) as the boundary of Lake Zhongshan, inside of which the high confidence segments were directly determined from the high RES reflectivity of the ice-water interface.
During the campaign two GNSS base stations were usually installed on the runway area to facilitate differential GNSS 235 data processing of the airborne GNSS positioning data. The aircraft GNSS receivers started recording at least 30 min before each flight to ensure pre-flight reference measurements (Cui et al., 2018;2020b). Thus, the survey line positioning error is 25 cm, which was given for flight lines surveyed using differential GNSS positioning (Cui et al., 2018). Furthermore, the RES sampling spacing is ~17 m. Combining these two error sources, we estimate that the uncertainty of the lake boundary uncertainty is ~17 m. The vertical uncertainty of the lake attributes are related to the 240 accuracy of the airborne GNSS positioning and the RES and gravity observations. An analysis of the ice thicknesses at the crossovers from different RES survey lines in PEL showed a vertical accuracy of ~10 to ~40 m (Cui et al., 2020b), which can be treated as the accuracy of the derived lake surface. Furthermore, the inferred lake depth has an accuracy of 50-70 m.
The lack of geophysical observations in PEL has resulted in uncertainties in modeling the subglacial lake and 245 channel system in this region (Mackie et al., 2020). The results of this study and the updated version of BedMachine (Morlighem et al., 2020) will provide important knowledge and data for the study of the mechanism of subglacial hydrodynamics that controls the water discharge from the interior of PEL to the Amery Ice Shelf and the Southern Ocean (Bell, 2008;Jamieson et al., 2016;Pattyn et al., 2016).
The same ELSA onboard data acquisition system of Snow Eagle 601 was tested at Lake Vostok on January 24, 250 2016 and collected data along a ~297 km U-shaped surveying line in the western part of the lake (Fig. 5a). The RES profile shows a clear distinction between brightness features in radargram from bedrock and the ~68 km long lake segment BB' (Fig. 5b). For consistency with the computation in Lake Zhongshan, the sediment area of higher reflectivity (left of point B in Fig. 5b) is excluded from the high confidence segment BB'. In accordance, we detected a maximum reflectivity increase of 23.8 dB and an average reflectivity increase of 8.6 dB from the bedrock to the lake 255 surface (Fig. 5c). Although the largest reflectivity increase measurements of all seven segments in Lake Zhongshan, 10.2 dB (maximum) and 5.1 dB (average) ( Table 1), are lower than those of the known subglacial lake environment in Lake Vostok, we notice that the reflectivity of the bedrock itself in the Lake Zhongshan area is also significantly https://doi.org/10.5194/tc-2021-332 Preprint. Discussion started: 3 December 2021 c Author(s) 2021. CC BY 4.0 License. lower than that of the Lake Vostok area (by ~7.0 dB). This may be attributed to the higher attenuation within the ice layers in the Lake Zhongshan region, as indicated in Figure 7 in Matsuoka et al. (2012). Therefore, in such a region 260 with highly attenuated RES signals through ice layers (Hills et al., 2020), the reflectivity contrast between water and bedrock may be hampered. Thus, the calculated average reflectivity increases in Lake Zhongshan, 5.1 dB (strongest segment) and 3.6 dB (average of all seven segments), are believed to be a safe indicator.

275
In summary, we developed an innovative and systemic method to build a set of evidence of the newly detected subglacial lake, Lake Zhongshan. We designed a targeted and intensified data acquisition strategy to collect both RES and gravity data in the PEL region using the ELSA system onboard the fixed wing aircraft, Snow Eagle 601, during the CHINARE expeditions from 2015 to 2019. Adaptive RES data analysis allowed us to detect the ice-water interface and the lake extent. We quantified the lake depth, volume, and other characteristics by gravity inversion modeling. 280 We report that Lake Zhongshan is located at 73°26'53"S, 80°30'39"E and is ~3603 m below the ice surface of PEL, East Antarctica. Its dimensions are 45,20028 m in length and 10,90028 m in width. The lake area is 3281.5 km 2 and thus, it is the only subglacial lake detected in PEL and by area it is the fifth largest in Antarctica. The deepest part of the lake bed is 25070 m from the lake surface and ~900 m below the sea level. The water volume of the two-basin lake is ~54 km 3 . These findings are important to understand surface ice dynamics and subglacial hydrodynamics in 285 PEL, as well as the stability of the Antarctic Ice Sheet. The results can also provide supporting data for potential drilling into the lake. In the future, the evolution of this subglacial lake will be explored with additional remote sensing data, continued airborne data collection, and eventual in situ or core observations. 290 Data availability. The RES data are currently included in appendix as graphyics. Readers can access the data at NSIDC (https://nsidc.org/data) when the paper is accepted.
Author contributions. RL and BS led the study, LL, RL, AZ, TF, LA, XC and TH carried out data analysis and editing, BS, XC and JG collected data, XC, BX, SL and LJ processed RES data, LA and LL performed gravity modeling, all authors were involved in writing and presentation. 295 Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. Logistic and technical support provided by the CHINARE program, the Polar Research Institute of China, and the Snow Eagle 601 operation team is appreciated. https://doi.org/10.5194/tc-2021-332 Preprint. Discussion started: 3 December 2021 c Author(s) 2021. CC BY 4.0 License.