An Unmanned Aerial Vehicle Sampling Platform for Atmospheric Water Vapor Isotopes in Polar Environments

Above polar ice sheets, atmospheric water vapor exchange occurs across the planetary boundary layer (PBL) and is an important mechanism in a number of processes that affect the surface mass balance of the ice sheets. Yet, this exchange 10 is not well understood, and has substantial implications for modeling and remote sensing of the polar hydrologic cycle. Efforts to characterize the exchange face substantial logistical challenges including the remoteness of ice sheet field camps, extreme weather conditions, low humidity and temperature that limits the effectiveness of instruments, and dangers associated with flying manned aircraft at low altitudes. Here, we present an Unmanned Aerial Vehicle (UAV) sampling platform for operation in extreme polar environments that is capable of sampling atmospheric water vapor for subsequent 15 measurement of water isotopes. This system was deployed to the East Greenland Ice-core Project (EastGRIP) camp in northeast Greenland during summer 2019. Four sampling flight missions were completed. With a suite of atmospheric measurements onboard the UAV (temperature, humidity, pressure, GPS) we determine the height of the PBL using on-line algorithms, allowing for strategic decision making by the pilot to sample water isotopes above and below the PBL. Water isotope data was measured by a Picarro 2130-i instrument using flasks of atmospheric air collected within the nose cone of 20 the UAV. The internal repeatability for δD and δO was 2.8 ‰ and 0.45 ‰, respectively, which we also compared to independent EastGRIP tower-isotope data. Based on these results, we demonstrate the efficacy of this new UAV-isotope platform and present improvements to be utilized in future polar field campaigns. The system is also designed to be readily adaptable to other fields of study, such as measurement of carbon cycle gases or remote sensing of ground conditions.


Introduction 25
The Greenland and Antarctic ice sheets interact with the atmosphere through continuous exchange of water vapor by condensation and sublimation, and through precipitation events (Fettweis et al. 2019). The planetary boundary layer (PBL, the lowest layer of the troposphere directly influenced by the surface) generally has a thickness of 10s to 100s of meters above the ice sheet, and exchanges water vapor with the free troposphere (FT) (Helmig et al. 2002). It is not clear how much water vapor is exchanged from surface sublimation flux, nor if the exchange ultimately results in a Bridging the two different scales of satellite remote sensing and in situ ground-based measurements is a challenging 65 necessity for understanding the hydrologic cycle. Most efforts and testing have occurred at lower-latitudes, far from the ice sheet. Franz and Röckmann (2005) developed a cryogenic sampler and protocol to collect stratospheric water vapor, from very small mixing ratios (<10 ppm) flown on a C-17 aircraft during flights between New Zealand and Antarctica. In 2007, Strong et al. was successful in using pre-evacuated 650 cc glass flasks to collect atmospheric water vapor samples in the field, then cryogenically extracting the water and reducing it to hydrogen (Friedman et al. 1953), followed by mass 70 spectrometer analysis. Vertical profiles were collected in approximately 300 m intervals using a light manned-aircraft with a ceiling of 2-3 km above ground level (AGL) in the desert southwest of the U.S. (Strong et al. 2007). As the engine of the aircraft was turned off during sampling in the Strong el al. study, obtaining airborne samples near the surface would be too dangerous.

75
There have been two recent measurement campaigns that utilized in-situ optical water vapor isotope instruments to constrain remote-sensing water isotope products. ISOWAT-II instrument over the Canary Islands to triangulate between ground-based Fourier transform infrared (FTIR) spectrometer measurement and space-based IASI (infrared atmospheric sounding interferometer) during the MUSICA campaign (MUlti-platform remote Sensing of Isotopologues for investigating the Cycle of Atmospheric water). A validation study estimated a 40‰ uncertainty of δD in the lower troposphere and 15‰ in the upper troposphere against the FTIR product. Uncertainty in IASI was estimated by Schneider et al. in 2015 to be 15‰ in the mid troposphere with a +30-70‰ 85 bias. Uncertainties of this magnitude are inadequate for constraining water vapor across the PBL and remain a target for improved methodologies.
We present results from a UAV pilot study at the East Greenland Ice-core Project (EastGRIP) site in northeast Greenland, occurring in summer 2019. We describe how customized UAVs can now be used to safely bridge satellite and ground-based 90 measurements, all while overcoming the challenging polar conditions to sample atmospheric air in the low-to-mid troposphere above the Greenland Ice Sheet. This is accomplished by designing an effective yet relatively inexpensive sampling platform with 3D-printed parts and accessible control devices on a commercially available fixed wing UAV that collects air samples aloft for analysis immediately following flight with ground-based instrumentation. We show that water vapor isotope measurements can be achieved with sufficient precision relative to the magnitude of the observed gradient 95 across the PBL, and comparable with independent measurements made at the EastGRIP 10m tower. We also demonstrate that algorithmic methods of evaluating clustering indices of real-time on board sensors to determine the altitude of the PBL, which can be used by the flight team to make informed sampling decisions mid-flight. We make recommendations for future field deployments to polar ice sheets and discuss the potential for how the observations can be used to improve the scientific understanding of varying fields of study. 100 2 Methods

Water Isotope Measurements
In this study, we made atmospheric water vapor measurements at the EastGRIP ice core field site in northeast Greenland (75.63°N, 35.99°W; 2,700 m above sea level). A cavity ring-down laser spectroscopy (CRDS) instrument, model L2130-i (Picarro Inc., Santa Clara, CA) was used in conjunction with a custom inlet to introduce both samples and standards with 105 equal treatment, described in more detail in Section 2.6. The standard water isotope data was analyzed on a continuous flow analysis (CFA) system adapted from Jones et al. (2017a). Results were validated against measurements made by the SNOWISO project (H2020 European Research Council Start Grant #759526), also using a Picarro L2130-i instrument (Section 2.3).

110
The data consist of measurements of hydrogen and oxygen isotopes in water vapor, where the ratio of heavy to light water isotopes in a sample is expressed in δ notation (Epstein et al. 1953, Mook 2000 relative to internationally recognized primary reference materials Vienna Standard Mean Ocean Water (VSMOW) and normalized to Standard Light Antarctic Precipitation (SLAP) in accordance with IAEA reference material (2017): where R is the isotopic ratio 18 O/ 16 O or D/H (i.e., 2 H/ 1 H). The δD and δ 18 O symbols refer to fractional deviations from VSMOW, normally expressed in parts per thousand (per mille or ‰). In practice, we maintain a suite of secondary reference waters that are rigorously calibrated to the primary reference materials (VSMOW and SLAP). Storage of our secondary reference waters is in accordance with methods described in IAEA Technical Note No, 43, (Newman et al. 2009).

EastGRIP Hydrological Cycle 120
The hydrological cycle on the Greenland ice sheet has several isotopic reservoirs and exchanges ( Figure 1). The dominant reservoir is the ice sheet, composed of ice, firn and snow with a relatively positive water isotope value compared to the overlying atmosphere (Steen-Larsen et al. 2011). At the ice sheet/atmosphere interface, both radiative (shortwave and longwave) and non-radiative (sensible and latent heat) energy fluxes occur, affecting the energy mass balance of the ice sheet. The summation of these processes leaves a diurnal imprint on the water isotopes in the upper few centimeters of the 125 firn (Ritter et al. 2016, Madsen et al. 2019, Hughes et al. 2021. Within the PBL, turbulent mixing occurs with a magnitude largely dependent on stratification and wind shear. Significantly stable stratification of the PBL (e.g. during polar nights) may serve in part as a preventative mechanism of vapor leaving the ice sheet (Berkelhammer et al. 2016). At a constantly varying height above the ice sheet (10s to 100s of meters in summer, lower in winter), a mixing zone between the surface and the PBL-free troposphere boundary allows for entrainment of water vapor from the free troposphere 135 into the PBL. This exchange is not well understood due to the inability thus far to make measurements across the full PBL (Boisvert et al. 2016). The inclusion of outside air parcels is mediated by synoptic changes in atmospheric general circulation (Schuenemann et al. 2009). Characterization of these synoptic scale changes have been shown to be important to large scale melt events, such as the 2012 event across the Greenland Ice Sheet where changes in atmospheric circulation resulted surface melt (Hanna et al. 2014). Due to the conservation of water isotopes through mixing, gradients in water isotopes across the 140 PBL-free troposphere-mixing zone may provide evidence of the amount of water vapor exchange between air parcels. As UAV methodologies improve, it will eventually be possible to provide constraints on net exchange of water vapor across the PBL-free troposphere interface.

EastGRIP tower measurements
During our UAV field campaign, simultaneous measurements of water isotopes were continuously taken at several heights 145 above the snow surface. The tower set-up used for these measurements was similar to the system described in Madsen et al. (2019). Four air intake inlets were installed at 0.5, 1.0, 2.0 and 7.1 meter height above the snow surface from which air was pumped to a Picarro L2140-i analyzer in a temperature-controlled tent ~15m away using an auxiliary pump.
In addition to documenting a diurnally varying water vapor isotope signal, the tower measurements have successfully been 150 used to observe a gradient in the isotopic concentration in the lowest part of the PBL (Ritter et al. 2016, Madsen et al. 2019.
This gradient has been used to argue that the exchange between the atmosphere and snow surface is driving the diurnal water isotope variations. Extending beyond tower heights will allow for the observation of entrainment processes and a better understanding of the formation of the ambient isotopic composition.

Fixed Wing UAV Flight System 155
While at the EGRIP camp in 2018 the team performed a proof of concept for airborne sampling and surface analysis using a small remote controlled sampling package and a multi-rotor UAV (DJI S-1000, DJI, Inc.) The system was able to obtain data and samples for analysis up to 400 meters AGL, but navigation and control was very problematic, due to proximity to the magnetic pole and batteries at low temperature limited flight times to less than 15 minutes. Knowing that sampling was possible and effective, we moved our attention to fixed wing platforms that fly longer, higher and are more stable to operate. 160 The S2 fixed-wing aircraft was the chosen platform for the 2019 campaign. The S2 is a modular, autonomous, aircraft designed by Black Swift Technologies, LLC (BST) for science missions, based on simple to operate electric propulsion aircraft with a modular payload. It includes a lightweight composite airframe design ( Figure 2). The S2 is capable of conducting fully autonomous flights in unimproved areas such as an ice sheet in part due to its pneumatic launch system. 165 The aircraft can adjust to changing wind conditions in real-time, ensuring a high degree of stability for predefined mapping or atmospheric sampling applications (Elston et al. 2015b). The aircraft can carry up to a 3.5 kg payload for up to 90 mins.
At arctic temperatures with the payload used in this study, we found 45 mins of flight time typical and apt for climbing 1600m and including needed sampling time. The broader list of technical specifications for the S2 are listed in Appendix C.
A typical flight day including sampling is found in Section 2.8. 170 SwiftPilot™ (Black Swift Technologies, Boulder, CO) is a miniaturized autopilot system developed specifically for UAV applications, allowing for remote operation and autonomous operation monitoring with capability for intervention, and was used in this study. Its modular CAN-bus architecture enables a large number of connectivity options, simplifying payload integration into the processing stream. Communication with the ground is enabled through the SwiftStation™ (BST), a portable tripod-mountable ground station (1.8 kg) that supports user-specific sensor payload integration, downlink, waypoint programming and digital terrain model custom inputs, and operation control. The standard configuration, used in this study, contains a 3dBi gain 900 MHz dipole as well as a GPS antenna. 180

Nose Cone Sampling Pod
The flask sampling apparatus is contained within the nose cone, and a schematic of the system is shown in Figure 3. The payload is suspended on four carbon fiber rods spaced 140mm x 80mm apart which slide into the frame of the main aircraft where a manufacturer-supplied baseplate secures it in place with two spring-loaded latches. We explored and tested the efficacy of holding water vapor within Teflon, Tedlar, and stainless-steel bags and we observed memory effects in all three 185 of those options. Glass was the only material found where sample carry over was minimal. As such, eight glass flasks (Precision Glassblowing, Denver, Colorado) are suspended with memory foam in a series of modeled and 3D printed nylon-12 plates (KODAK Nylon 12). Due to 12G launching force from the pneumatic launching process, we found foam and the elastic properties of nylon-12 to be critical for flask safety. The printing was done on a XYZprinting da Vinci Super and

Measurement Scheme
Similar sample pod control systems were used for both airborne and ground sampling. For sampling during flight, the onboard microcontroller (Adafruit Feather M0) works through the BST SwiftCore™ flight system to communicate to the ground station. Payload control is managed by a laptop with Linux (Ubuntu 18.04.2) connected over WiFi to the ground station. The microcontroller receives and manages commands to toggle valves and enable pumping. Environmental sensing 210 is also fed into the BST SwiftCore™ and down to the ground station. The temperature and humidity is determined by an E+E Elektronik EE03 sensor (±0.3°C and ±3%RH), and the pressure is determined by a high resolution (±1.5mbar) MEMS sensor (TE Connectivity MS5611). Both sensors are included as part of the forward pointing package to assist in autopilot flight on the right wing of the aircraft.

215
In addition to measurements of samples taken during flights, a small (2 meter) sampling tower was used for flask sampling to provide an additional near surface data point and also allow an intercomparison with tower measurements of water vapor isotopes at EastGRIP. On the ground, a second microcontroller was connected to the sample pod with a USB cable. Its tasks included controlled functions 1) flushing dry air through flasks prior to flight, 2) sample acquisition from the 2 meter tower, and 3) computer controlled release of samples for isotopic analysis. Flasks from both flights and ground sampling are 220 introduced to an L2130-i Picarro instrument for isotopic analysis by opening a single port on the flask. Before air sample ports are opened, dry air is plumbed into a spare valve at the back of the manifold to push out atmospheric air left in the manifold. In this manner the Picarro analyzer is pulling the sample air from the dead end of the flask, reducing the pressure slowly over time. Air samples are pulled from the flask into the instrument via the common port of the manifold at 30 sccm for approximately 12 minutes. Pressure within the analyzer cavity is carefully controlled at 50 Torr by the instrument with 225 high speed PID controlled valves on both ends of the cavity. As water vapor is introduced to the CRDS cavity, isotopic mixing with the previous dry air parcels can affect the instrument's response to new samples. To address this, the first 3 minutes of observation for any one sample is cropped from averaging. Additionally, to address any issues associated with any reductions in flask pressure near the end, the last 3 minutes are also cropped. These timings were empirically derived from consistent plateaus of both isotopes and water concentrations between the beginning and ending tails. Cropping in this 230 way also allows a mixing ratio/specific humidity to be determined for calibration. Values for any one sample are determined from the average over approximately 6 minutes. For a systematic diagram of the drone and ground sampling, see Appendix A.
The methods insured equal treatment of samples collected in-flight or on the ground. This served two purposes, 1) to 235 establish the isotopic bottom end-member of the vertical profile. and 2) to enable the comparison of the sample pod measurements with the established in-situ tower measurements of water vapor concentration and isotopes (Picarro L-2140-i), taken at the same time within a distance of 10 m.

Water Vapor Isotope Measurements and Calibration
Systems have been developed by numerous groups to calibrate Picarro CRDS instruments used in continuous flow 240 applications (Gkinis et al. 2011, Jones et al. 2017a, and each represents an evolution in design and performance. Due to the proven success with multiple measurement campaigns completed on ice cores with the calibration setup described in Jones et al. (2017a), we used the same principles in this setup for the calibration of the system in the field.
It meets the ideal criteria for a calibration system as described in Bailey et al. (2015), that includes a) enabling the introduction of low volume mixing ratios for calibration, b) mitigating standard drift, and c) utilizing multiple water 245 standards in the calibration scheme. The system schematic is shown in Figure 4.

265
Raw values from the CRDS system are corrected in post processing and tied to known values of isotopic water standards and corrected for the instrument's response to humidity. This is required because our atmospheric water vapor samples (typically <5,000 ppm H 2 O) are outside of the standard operating range of the Picarro L-2130i, which is optimized for the analysis of liquid water samples (10,000 to 25,000 ppm H 2 O). Counting statistics for CRDS instruments are heavily dependent on sufficient concentration of gas species warranting calibration across a range of humidities and isotope standards. The 270 isotopic water standards and their uncertainties are given in Table 1.  This calibration procedure was done several times throughout the 2019 field season to capture long-term instrument noise in response to humidity. Atmospheric samples were calibrated to the set of humidity measurements closest in time, ranging 290 from as long as 7 days apart but typically 1-3 days throughout the season. Steen- Larsen et al. (2013) indicates that correctable linear drift may occur local in time to the measurement period due to strong diurnal temperature changes around the instrument. Because humidity calibrations were not regular about each measurement at the time scale of diurnal temperature change, the correction was not performed in this study. Future campaigns will include a higher calibration density to account for this. 295

Uncertainty in Sampling and Intercomparison with On-Site Water Vapor Tower
Outside of CRDS instrument performance, the UAV sampling system itself introduces sources of error. This uncertainty is associated with acquisition and transport of the sample water vapor as well as environmental change during the flight period.
To understand the uncertainty in captured water vapor during the 2019 season, two different flask pod intercomparisons were performed in conjunction with the separate 2-meter tower-isotope setup detailed in Section 2.3. For the intercomparison, 300 each of the six flasks from three different sample pods was flushed with air from 2-meter altitude for 5 minutes. As a total of eighteen flask measurements correspond with an hour and half of sampling, this test is sensitive to changes in atmospheric water vapor isotopic composition. A more appropriate test would be to produce standardized water vapor as described in Section 2.5 and sample from that stream. This is challenging because the most accurate test would be to produce water vapor at a rate that can match the 5 LPM sampling throughput of the pump, which is currently unachievable due to limited amounts 305 of water standards. Though sampling was performed over this longer period of time without standard water vapor, the highest 1σ standard deviations of any one pod was of 0.45 ‰ in δ 18 O and 2.80 ‰ in δD. These values can be seen as the pessimistic view of uncertainty due to the non-ideal sampling situation, but are reasonable given that previous uncertainty estimates on in-situ water vapor isotope measurements range from 0.14 ‰ in δ 18 O and 0.85 ‰ in δD (Steen-Larsen et al. A comparison of UAV and tower deuterium excess (dxs) data is shown in Figure 5. Deuterium excess is defined by Dansgaard (1964) as dxs = δD − 8δ 18 O. The dxs is a more sensitive intercomparison metric than δ 18 O and δD and will more clearly show discrepancies between different measurement schemes. An intercomparison was done at four different times: 1) during 2 meter sampling during two different flights and 2) during two different pod intercomparison measurements at 2 315 meter. There is general agreement for dxs across the two platforms with a slightly more positive value for the UAV-isotope system. The positive relation is seen in both δD and δ 18 O implying that the positive bias is due to an interplay of both measurements. Figures of separated δD and δ 18 O can be found in Appendix A.

Boundary Layer Prediction
During the 2019 field campaign, we used environmental measurements (pressure, potential temperature, specific humidity) 325 taken in real-time during each flight to evaluate Euclidean distance in the measurement domain to infer where the PBL/free troposphere transition occurs in the spatial domain (Appendix D). The results were used by the pilot to make in-flight decisions about sampling altitudes for isotopic analysis. After the 2019 field campaign, we explored additional PBL identification algorithms. The PBL and free troposphere are largely decoupled, allowing for cluster density evaluation to determine the PBL height (Krawiec-Thayer 2018). As the PBL structure varies in shape and magnitude for any one 330 observational parameter, other methods such as gradient interpretation of single environmental variables are less useful (Krawiec-Thayer 2018). The most promising algorithm, the Calinski-Harabasz criterion index (CHCI), is explained in Appendix D. The global maximum of this index is assumed to be the height of the PBL. The Calinski-Harabasz criterion index will be utilized in future field campaigns to detect the PBL in real-time during flight in addition to user judgment. An example of the CHCI for PBL height determination is shown in Figure 6. 335 Figure 6: The Calinski-Harabasz criterion index (CHCI) applied to the sampling flight on June 12th with a-priori assumption of K=2. Groups are assumed to represent the free troposphere and PBL, though more structure may exist. The boundary at ~325 m using CHCI is not in agreement with that determined using the Euclidian distance in Figure 7, which shows a likely boundary at ~275 m. This value was chosen by the operator of the flight and shown as the dashed green line. Figure 10 shows an example of a 340 failure in Euclidian distance to predict the boundary layer. The use of CHCI improves the PBL prediction algorithm, as determined in this study.

Typical analysis day and sample acquisition
Before flight is considered, the local weather is evaluated to determine the potential for mission success. To prevent potential icing, a nearby ceilometer (Vaisala Ceilometer CL31, Vaisala, Boulder, Colorado) present at the EastGRIP camp was used to 345 safely determine that cloud cover was significantly higher than the highest flight altitude in the flight plan. Flights were not performed during precipitation events. Acceptable wind speeds were considered less than 10 m/s, two-thirds of maximum wind operation of 15 m/s for the Black Swift S2 aircraft.
For any given analysis flight, a sequence of steps are completed to ensure quality control: 1) Calibration of the water isotope 350 measurement system (Section 2.6), 2) On going isotopic measurements at a 2-m tower during the flight (Section 2.7), 3) Identification of the PBL during flight using real-time temperature and R/H from the aircraft (Section 2.8), 4) Atmospheric sample acquisition during flight, and 5) Isotope measurement following the flight, in a heated field tent (Section 2.6).
A calibration of the Picarro L-2130i is performed close to the time of flight. Before a flight, both ground-based and UAV-355 based glass flasks are flushed with dry air (75 ppm water vapor) for 10 min. Before launch (time permitting), an extra 2 meter measurement is taken with the ground sampling system detailed in Section 2.5.2. After launch, the pilot ascends at an autopilot-controlled rate of 2 m/s in a circular pattern (a 68 m diameter orbital). The ascension rate can be affected by local wind speeds requiring a slower vertical climb than the UAV is otherwise capable of. While a faster ascension is possible, a slower climb also minimizes hysteresis for the atmospheric sensors onboard the UAV. At the top of the climb, the aircraft 360 automatically enters a holding orbital pattern at constant altitude while the operator assesses the real-time algorithmic determination of the PBL. The operator then inputs the altitude of the sampling locations for water isotopes above and below the PBL.
The UAV then descends to the first/highest sampling altitude. At each sampling altitude, the pilot initiates flask sampling. 365 The sample procedure can be broken into three steps: 1) holding altitude, 2) flushing, and 3) equilibration. When the UAV reaches the first sampling altitude, the UAV will maintain altitude for approximately one minute to eliminate hysteresis of the environmental sensors. The diaphragm pump is then turned on and each port on the flask is opened for a three-minute flush of ambient air to address memory effects on the interior glass surfaces. Then, the pump is turned off in order for the flask to equilibrate to ambient pressure for 10 seconds. Finally, the valves are closed, and the process is repeated for a second 370 flask, providing paired measurements at each altitude. Paired sampling was motivated primarily by the inability to test the low temperatures, the 12G forces exerted on the flasks during launch, and inflight vibration forces in a "benchtop" setting.
The nose cone sampling pod holds 8 flasks, allowing for paired measurements at four altitudes. However, due to battery limits on site, the payload was generally flown with 6 flasks (3 pairs). The aircraft is then directed to land. Both the UAV atmospheric samples and ground-based samples (from 2 meter height) are then analyzed on the water isotope measurement 375 system and calibrated to the most recent system calibration (Section 2.5).

Retrieval of Water Vapor Isotopic Composition about the PBL
Though CRDS measurement of water vapor isotopes by aircraft is not new (Section 1), its capture and retrieval by UAV for later measurement is novel. Arctic environments present major logistical challenges for fieldwork. The remoteness of field 380 camps, such as EastGRIP, makes logistics challenging and limits the amount of field personnel. The potential for extreme weather, cold temperatures, blowing snow, and safety are all significant factors that limit scientific outcomes. For these reasons, even the most careful planning will still result in some unforeseen challenges. During our field campaign, we realized that we had to improve system sampling turnover time to produce more flights per day, that hysteresis in the environmental sensor could produce artifacts in PBL detection, and that our 2-3 person field crew was inadequate to have 385 good diurnal sampling coverage since all people slept during the same hours. A larger team would have provided an option for day and night shifts as there were 24 hours of sunlight during the field campaign.
Despite unforeseen challenges, we achieved a total of four sample-taking flights from June 12th to June 26th, 2019. An example of environmental sensor data for June 12th is shown in Figure 7. We found varying amounts of structure in isotope 390 space across all four flights (Figure 8). Large transitions between water vapor isotope surface measurements at 2 meter and values above and below the PBL/free troposphere (FT) transition are apparent in the June 12 th flights (Figure 8). The other flights in contrast had little variability, suggesting that the PBL was unstable (i.e. well mixed). Berkelhammer et al. (2016) suggested that summertime nights at Summit, Greenland would present the conditions for stable stratification of the atmosphere, but that this claim was unprovable using towers alone. In 2022, we will use an improved UAV-system setup to 395 generate a comprehensive diurnal data set spanning many weeks' worth of time.  a and b)

Hysteresis and Calinski-Harabasz Criterion Index (CHCI) and PBL Detection 410
The CHCI was calculated post-flight for comparison with 1) the self-similarity of Euclidean distance (used during the 2019 field campaign, but later updated to the CHCI approach) and 2) operator determination of the PBL. The results are shown in Appendix A. The CHCI had a direct match with Euclidean distance for half of the flights. In the other half, the CHCI predicted altitudes significantly higher than the other determinations. The results of our comparison reveal that our original PBL-detection algorithm using Euclidian distance needs improvement (Figure 9). Specifically, we have determined that Euclidean distance can under or overestimate the height of the PBL due to sensor (temperature and humidity) hysteresis.
This hysteresis exceeded the stated manufacturer response time for the atmospheric temperatures we encountered, discussed in Appendix B. The hysteresis is almost entirely a result of the rate of ascent during flight. Before a flight, the UAV is static at ground level, thus temperature and humidity measurements will be stable, varying only slightly with small changes in surface conditions. The energetic pneumatically-driven launch of the aircraft (a 12 G force) results in a rapid increase in 420 altitude that can introduce a bias into the sensor output due largely to the thermal mass of the sensor and slow response to rapidly changing conditions. A similar effect occurs anytime the rate of ascent is not constant, such as when the UAV transitions between different orbitals (i.e. a sampling orbital and landing orbital).
A case study in Figure 9 illustrates a shift in orbitals from the June 21st mission. The operator moved from the initial launch 425 orbital to a lower altitude to begin an ascension profile. During the transition between the two orbitals, the aircraft moved from about 110m to 60m in altitude in ~1 minute. During the transition and immediately during the ascent, multiple temperature and humidity values were generated for the same altitude creating a region of varying hysteresis effects that can bias PBL prediction by Euclidean distance, ultimately causing the operator to misidentify the altitude of the PBL. More concisely, the algorithm detected this data anomaly as atmospheric structure, when in fact it was due to hysteresis. While 430 removing this skewed data could be an easy fix, the stabilization of temperature and humidity to that new starting altitude biases the beginning of the climb just as it does at the surface before launch. The hysteresis effect is also noticeable in the CHCI (Figure 10, green circles). Relaxing the a priori assumption of a single PBL that separates the surface atmosphere from the free troposphere, additional transition regions can be identified. As 440 CHCI uses Euclidean distance to establish variances, it is also subject to potentially poor predictions in situations of significant hysteresis. However, its ability to establish regions of similarity, such as the case of the transition region between launch orbital and the ascension orbital during the June 21th mission provides an objective method of informing the operator of potential false positives for the boundary layer altitude. In this specific case, three of the top five PBL altitudes predicted by the Euclidean distance algorithm can be flagged as incorrect. However, even with sensor hysteresis, we determine the 445 CHCI to be an effective tool to assist in fast mid-flight evaluation of the boundary layer by the drone operator. Figure 10: Specific humidity over the ascension for the June 21st flight partitioned into groups by the CHCI with the a-priori assumption of K relaxed from 2 to 4. The region of transition the operator took post launch between ~110m and ~60m is clearly evident as a separate group (green circles). In cases where artificial structure exists due to sampling patterns, CHCI may assist the 450 operator by flagging those areas.
Overall, there are two options for overcoming the effects of hysteresis: 1) better sensors and 2) changes to flight mission plans. We have identified the Vaisala RSS-421 sonde sensor to meet the first requirement. The RSS-421 includes a low thermal mass fine-wire thermocouple and heated humidity sensor with bakeout unit, which will allow for faster response in arctic conditions. This sensor has already shown to be capable of producing accurate temperatures in challenging UAV fixed-wing missions (Frew et al. 2020). For flight planning, relocating launch sites to be as close to the ascension orbital as possible will reduce hysteresis during horizontal transitions between orbitals. The ascension rate can also be slowed to less than 2 m/s allowing the maximum time for sensors to equilibrate with the surrounding atmospheric conditions. The tradeoff is that this may require reducing the maximum flight altitude to conserve battery life and reduce the bank angle. A sharp bank angle decreases the lift coefficient (Williamson 1979), and a higher angle of attack is needed to maintain ascension rate 460 in tailwind situations (Blakelock 1991). When the pitch angle needed is too high and outside the flight envelope, the Black Swift Technologies autopilot will slow ascension to protect the aircraft. It is assumed that variability in temperature, pressure, and humidity is small in the x and y plane, allowing for a large increase in orbital diameter to reduce bank angle significantly.

Conclusions and Outlook 465
We have presented a UAV-isotope sampling platform and methodology capable of measuring atmospheric water vapor and its stable isotopes within the planetary boundary layer (PBL) and lower troposphere in a polar environment. We utilize a fixed-wing UAV (Black Swift Technologies) with flight times in excess of 45 minutes with the capability to reach 1,600m AGL. Multiple nose cones allow for collection of air in 8 glass flasks, enclosed within a 3D printed support structure that critically withstands 12Gs of force during takeoff. In this study, the total system is used to sample above and below an 470 algorithmically-detected PBL, resulting in the first measurements of atmospheric water isotopes above and below the PBL on the high-altitude Greenland Ice Sheet.
Across four sample-taking missions at the EGRIP ice core site in 2019, we observed significant variation in water isotopes on either side of the PBL; the variability exceeded our conservative precision estimates of 2.8‰ in δD and 0.45‰ in δ18O. 475 These results form the basis for future campaigns to collect high-temporal density measurements (flights every 4-6 hours across many weeks) at key missing scales that will improve ice-to-atmosphere modeling and mixing processes, ice sheet mass balance, satellite detection algorithms, moisture tracking, ice core science, and modeling the hydrologic cycle in general.

480
A field campaign for return to EastGRIP is scheduled for summer 2022. Future improvements to the UAV-isotope system will be primarily focused on logistical improvements that increase the number and frequency of flights. Additional flight crew will be available for nighttime flight missions. To ensure a balanced diurnal flight schedule over weeks of time, with the goal of one flight every 4-6 hours, a precessing schedule of calibration times will be used. Each calibration will be done every 2-4 days, lasting 12 hours, starting at different times of day. This ensures that we do not consistently lose the ability 485 for UAV sampling at the same time for every calibration, e.g. from 12pm-12am. The combination of these improvements will allow the potential maximum number of flights per day to increase from two to as many as six, while balancing the timing of calibration. In flight, we will carefully regulate the rate of ascent and include better performing temperature and humidity sensors with minimal time constants, all of which will reduce hysteresis for PBL detection. We plan to leverage an existing anemometer used by the autopilot in order to assist in the correction as well as produce an additional 2D wind speed 490 for the flight. Additional improvements will include a lighter pump and manifold system that should allow greater flight time. Beyond Greenland, this platform is readily adaptable to other scientific disciplines, and will be used in an upcoming permafrost project to measure atmospheric methane emissions and soil moisture content in Alaska.       A supplier listed response time for temperature measurement of the ee03 sensor was not available within the temperature ranges measured within the study and assumed to be negligible.

Appendix C: The S2 Drone
Scientific missions the S2 has flown prior to this study include mapping soil moisture with a radiometer (Dai et al. 2016), a 535 calibration mission including a 12-band multispectral camera system (Wang et al. 2016), measuring snow-water equivalent with a radiometer (Yueh et al. 2018), and a volcano sampling mission that involves difficult operations into the plume of an active volcano (Wardell et al. 2017). The S2 is currently in use by the National Oceanic and Atmospheric Administration (NOAA) for wildfire applications (Gao et al. 2017) and it has flown in various challenging environments including at high altitude during atmospheric sampling campaigns in the San Luis Valley in Colorado (de Boer et al. 2018). The S2 is 540 designed for operations at altitudes up to 6,000 m AMSL in support of The National Aeronautics and Space Administration (NASA) science missions .
The S2 utilizes the SwiftCore™ Flight Management System for avionics control, communication, and command, designed by BST. It comprises the SwiftPilot™, SwiftStation™, and SwiftTab™ user interface, along with support electronics. 545 SwiftTab™ runs on Android devices like smartphones or tablets. Flight plans 1) can be uploaded, created, and modified before and during flight, 2) can use georeferenced data points for systematic surveying including pre-defined banking and spirals, and 3) are fully autonomous from launch to landing. Immediate preliminary analysis and decision making is supported via real-time telemetry and control capabilities. While Euclidean distance is more robust than individual gradient analysis (Krawiec-Thayer 2018), the technique still returns 565 multiple candidates for the PBL height. Instead, indexing methods can provide a deterministic global maximum of centroid partitions associated with the dataset. For the Calinski-Harabasz criterion index, the centroid is determined with a nonhierarchical k-means method. k-means is a data-partitioning algorithm that determines groupings of k amount of centroid clusters of n total observations converging to a maximum criterion value or index between centroids. k is determined a priori to be 2 corresponding to the assumed present atmospheric regions, the PBL and free troposphere. The Calinski-570 Harabasz criterion index has been used successfully with k-means methods in previous remote sensing and weather balloon studies (Toledo et al. 2014, Caicedo et al. 2017).
The Calinski-Harabasz index is the ratio of variance within one centroid and the variance between origin locations of all other centroids. Let m i as the centroid of cluster i containing n i data points, and c be an origin point for the data set. The 575 variance within one cluster is defined below in Equation 3: The expression for variance between clusters is defined as The ratio of variances, the Calinski-Harabasz index, then follows as 580 The centroid pair with the highest index is then the most significant group of partitions and the height that corresponds with the boundary of the two groups is assumed to be the upper layer of the PBL. An example of this method is shown in Fig. 6 and 10.

Data Availability 585
A data package upload has been initiated with Arctic Data Center, which is committed to providing citable datasets to facilitate reproducible science. Each DOI issued by the Arctic Data Center is intended to represent a unique, immutable version of a data package. We will finalize the data package during the review process with AMT, based in part upon referee feedback for the manuscript.

Acknowledgments 590
The authors are grateful for the funding provided by the National Science Foundation Award 1833165, "Closing the Water Vapor Exchange Budget Between the Ice Sheets and Free Atmosphere", managed by Jennifer Mercer. We wish to thank Dorthe Dahl-Jensen, University of Copenhagen and the EastGRIP international team for their support of the fieldwork on the Greenland ice sheet; EastGRIP is directed and organized by the Centre for Ice and Climate at the Niels Bohr Institute,

Author Contributions
HCSL, BHV, and TRJ developed the initial idea, rationale, and experimental setup. KSR, TRJ, and BHV prepared the original draft and all authors contributed to the review and editing of this paper. Design of the UAV sample pod was a product of BHV, KSR and JE. Flights and field water isotope analysis were done by KSR, VM, BHV, TRJ, WS, and AH. 610 Boundary layer prediction algorithms were developed by TRJ and KSR. Figures were prepared by KSR, BHV, and TRJ. For comparison to UAV flask samples, concurrent water isotope tower measurements were provided by HCSL, SW, and AKF.
Insights for modeling and error analysis were provided by HCSL, SW, and AKF.