Sea ice melt pond bathymetry reconstructed from aerial photographs using photogrammetry: a new method applied to MOSAiC data

Melt ponds are a core component of the summer sea ice system in the Arctic, increasing the uptake of solar energy and impacting the ice-associated ecosystem. They were thus one of the key topics during the 1-year drift campaign Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) in the Transpolar Drift 2019/2020. Pond depth is a dominating factor in describing the surface meltwater volume; it is necessary to estimate budgets and used in model parameterization to simulate pond coverage evolution. However, observational data on pond depth are spatially and temporally strongly limited to a few in situ measurements. Pond bathymetry, which is pond depth spatially fully resolved, remains unexplored. Here, we present a newly developed method to derive pond bathymetry from aerial images. We determine it from a photogrammetric multi-view reconstruction of the summer ice surface topography. Based on images recorded on dedicated grid flights and facilitated assumptions, we were able to obtain pond depth with a mean deviation of 3.5 cm compared to manual in situ observations. The method is independent of pond color and sky conditions, which is an advantage over recently developed radiometric airborne retrieval methods. It can furthermore be implemented in any typical photogrammetry workflow. We present the retrieval algorithm, including requirements for the data recording and survey planning, and a correction method for refraction at the air–pond interface. In addition, we show how the retrieved surface topography model synergizes with the initial image data to retrieve the water level of individual ponds from the visually determined pond margins.

We use the method to give a profound overview of the pond coverage on the MOSAiC floe, on which we found unexpected steady pond coverage and volume. We were able to derive individual pond properties of more than 1600 ponds on the floe, including their size, bathymetry, volume, surface elevation above sea level, and temporal evolution. We present a scaling factor for single in situ depth measurements, discuss the representativeness of in situ pond measurements and the importance of such high-resolution data for new satellite retrievals, and show indications for non-rigid pond bottoms. The study points out the great potential to derive geometric properties of the summer sea ice surface emerging from the increasingly available visual image data recorded from uncrewed aerial vehicles (UAVs) or aircraft, allowing for an integrated understanding and improved formulation of the thermodynamic and hydrological pond system in models.

1 Introduction

Melt ponds are a key driver of the summer energy budget on the sea ice surface. Their tremendous impact on the surface albedo and related self-reinforcing feedbacks lead to increased uptake of solar radiation (Fetterer and Untersteiner, 1998). However, the effects of melt ponds used to be parameterized rather simplistically in global climate models due to limited reference data, coarse resolution, and computing power (e.g., Pedersen et al., 2009a). Observational reference data that allow an integrated understanding of the thermodynamic and hydrological pond system are still rare (Wright and Polashenski, 2018). In particular, most melt pond depth observations used so far for model developments have been collected manually on the ice during comprehensive field campaigns, e.g., the Seasonal Sea Ice Monitoring and Modelling Site (SIMMS) (Morassutti and Ledrew, 1996) and the Surface Heat Budget of the Arctic Ocean campaign (SHEBA) (Perovich et al., 2003). Morassutti and Ledrew (1996) report greater pond depths on multi-year ice (MYI) (27.4 ± 12.6 cm) and landfast ice (LFI) (31.0 ± 19.2 cm) in comparison to first-year ice (FYI) (13.0 ± 8.0 cm). High standard deviations in the depth measurements were found between ice types and across spatial scales due to the inconsistent morphological nature of the different ice types. During the 1-year Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) from 2019 to 2020, Webster et al. (2022) found average depths of 22 ± 13 cm with manual measurements along transect lines. However, the actual pond bathymetry, which we define here as the pond depth profile in all directions and which therefore also yields the actual average pond depth, remains largely undiscussed in the literature.

Pond depth as a bulk property is used as a parameter in melt pond schemes in the Community Climate System Model version 4 (CCSM4) (Holland et al., 2012) and the implemented Los Alamos sea ice model CICE (Flocco et al., 2012). Holland et al. (2012) directly relate the use of different optical property parameterizations of the sea ice surface to pond depth and retrieve pond fraction from the available meltwater volume by formulating a linear relationship between pond fraction and pond depth from the SHEBA data, which gives them the filled pond volume. Pedersen et al. (2009b) developed a summer sea ice albedo scheme for the ECHAM5 general circulation model in which they derive pond fraction from pond depths given by surface melt rates. The link between fraction and depth in their scheme was developed based on a small-scale pond model by Lüthje et al. (2006).

A broader, more solid database of pond depths is still missing yet crucial for parameterizations in sea ice models building upon a deeper understanding of pond evolution and interactions. Promising new methods like the study by König et al. (2020) reveal that high-resolution optical remote sensing of the full pond bathymetry is possible on larger scales. They used the increased absorbance of radiation in liquid water at a wavelength of 720 nm to determine the thickness of the liquid water column independently of the pond bottom appearance from hyperspectral data. With this passive radiometric method, a spacious area could be covered by high-resolution optical data in the respective spectral band. However, the spectral method is restricted to observations under clear-sky conditions and is therefore still limited in application. Tilling et al. (2020) and Farrell et al. (2020) evaluated photon backscatter signals measured by ICESat-2 over sea ice and developed the UMD-MPA algorithm to derive the depth of particularly large (width >20 m) and deep ponds along the ground tracks of the satellite beams. Ongoing development of the ICESat-2 retrieval algorithms (DDA-bifurcate-seaice, in Herzfeld et al., 2023, minimum pond width of 7.5 to 15 m depending on the ice topography) and comprehensive data analysis presented in Buckley et al. (2023) highlight the ability to retrieve Arctic-wide pond depth data from satellites under cloud-free conditions. Another active technique to determine shallow-water bathymetry on a large scale is using airborne laser scanner (ALS) systems with water-penetrating wavelengths in the green spectrum as airborne laser bathymetry (ALB) systems. To our knowledge, such an ALB system over sea ice was only deployed for the first time in 2022 as part of the ICESat-2 validation.

Aerial RGB imaging platforms are numerously available and have been deployed in the Arctic for decades. They are already widely used to retrieve properties of bare surfaces in all different fields of geodetic studies. If the flight pattern is suitable, photogrammetric multi-view reconstruction can derive digital elevation models (DEMs) from aerial images. Sufficient forward and lateral overlap between images (about 80 % and 60 %, respectively) results in ground points being recorded from more than 15 different azimuth and elevation angles. This is used to achieve a triangulation-based reconstruction with a monocular camera system. A few studies could already show that the reconstruction method can be applied in mix-phased areas to retrieve the bottom topography of shallow riverbeds (Westaway et al., 2001) and coral reefs (Casella et al., 2017), as well as in laboratory seabed studies (González-Vera et al., 2020). From the nature of the method, the underlying ice or seafloor surface and some structure must be visible to be reconstructed. This limits the method to clear waters and shallow depths. Furthermore, appropriate correction methods for the light refraction at the water–air interface are needed. Although melt ponds probably closely conform to these requirements, a detailed method for deriving pond depth from aerial photographs has not yet been developed.

Only one experimental study on the photogrammetric derivation of pond depth from aerial images was carried out above sea ice before by Divine et al. (2016). They used a complex stereo-vision camera system on a helicopter to detect the sea ice surface morphology and melt pond depths north of Svalbard in 2012. Comparing their photogrammetrically derived ice freeboard results with terrestrial laser scanner data, they retrieved high agreement (bias of 0.03 m) and low deviation (rms of 0.04 m). Remarkably, they also retrieved the same accuracy when comparing melt pond depths with manually measured in situ data, although no physical correction was applied, which considers the different optical properties of water and air that make subsurface areas appear shallower. Potentially, their measured depths (<0.3 m) were too small to detect these effects. In contrast, Casella et al. (2017) measured significantly greater depths down to 1.8 m in coral reefs, also neglecting differences in optical properties, arguing that almost nadir measurements do not require corrections. Other studies on shallow-water bathymetry set up entirely new sets of equations to correct for light bending at the air–water interface directly in the photogrammetric reconstruction, requiring extensive efforts to solve the complex sets of equations in a reasonable time.

Given that photogrammetric methods for surface reconstruction from aerial images are well developed, reaching a user-friendly status, and the increasing availability of aerial images (also from drones) collected from monocular camera systems, we demonstrate here how corrections can be incorporated into a feasible workflow based on a commonly used photogrammetry suite to reconstruct entire melt pond bathymetry. Subsequently, we explore with the newly developed method how pond bathymetry evolved on the MOSAiC floe of leg 4 and how representative transect lines described the pond evolution on the entire floe, and we discuss possible upscaling factors for in situ point measurements.

2 Overview of data and ice conditions

2.1 Images and field data used in the method development from PASCAL 2017

We developed and initially tested the method on data from RV Polarstern Cruise PS106/1 (PASCAL) that took place in June 2017 (Macke and Flores, 2018) and subsequently confirmed it for deeper ponds (>0.4 m) with data from the MOSAiC expedition. During PASCAL, we collected aerial images (Fuchs and Birnbaum, 2021) together with coordinated manual ground-truth measurements. We measured in situ pond depths as a reference for the development of remote-sensing-based pond depth retrievals (including König and Oppelt, 2020; König et al., 2020) and as part of a geodetic survey of ponds, including their depth and bottom ice thickness. Coordinated data to study the photogrammetrical pond depth retrieval were available from 2 days: 10 and 14 June 2017. On both days, we collected RGB images above the in situ measurement area in a flight pattern that allows for the photogrammetrical retrieval of surface topography. In addition to a high overlap of images in the forward direction, which was present in nearly all images captured during the campaign, a crucial lateral offset of the acquisition positions was achieved on both days. The ice floe to which RV Polarstern was anchored during the campaign was located north of Svalbard at 81°50^′ N and 10°20^′ E. Given the location, we used zone 32N in the Universal Transverse Mercator system (UTM32N, EPSG:32632) as a projected coordinate system for all geospatial data evaluation.

Figure 1 Overview orthomosaics of the PS106/1 (PASCAL) study area north of Svalbard at 81°50^′ N and 10°20^′ E in June 2017. Marked ponds were examined in this study. White insets show the black reference targets enlarged. The known size of the targets was used for a scale check of the orthomosaics. The larger inset in (b) shows pond no. 2 used only in the 14 June 2017 evaluation. White arrows connected to the insets indicate their position at the study site, and meters give the distance from the edge of the orthomosaics. The maps are projected in UTM32N.

The study area was located approximately 1 km behind the stern of RV Polarstern in a younger, first- and second-year ice region that was, before our visit in June 2017, subject to strong deformation (König et al., 2020). Several depressions had been formed between rafted ice floes along a ridge and were flooded with seawater. At the time of our arrival, the PASCAL floe was already subject to melting, but visible melt ponds on level ice had not yet formed. The studied depressions were thus most likely initially formed by flooding and provided, at the time of our stay at the floe, a sink for incoming solar radiation, which then led to a catalyzed melt pond formation in the particular region. Due to the diverse appearance of the underlying ice, from bright blue to almost black pond bottoms, this was an ideal study area (Fig. 1).

Images on 10 June (Fig. 1a) were acquired from the RV Polarstern helicopter D-HARK with the implemented CANON EOS 1D Mark III 14 mm lens nadir camera system during measurement flight 20170610-2, which took advantage of the thoroughly clear-sky conditions and aimed at an upscaling of ground measurements. Therefore, several flight legs were flown at different flight levels (60 to 3000 m) above the study area. For this method development study, relying on high-resolution data, all images were used that have been captured below 200 m flight altitude and therefore provide a ground sampling distance (GSD) more precise than 10 cm, which is at least 1 order of magnitude smaller than typical pond extents at the study site.

Weather conditions on 14 June (Fig. 1b) were the opposite. The entire sky was covered by a stratiform cloud cover, as is usual in central Arctic summers (Cotton et al., 2011) when the average cloud coverage reaches its maximum of about 70% (e.g., Wang and Key, 2005). No solar disk was visible, meaning incident light can largely be assumed to be diffuse at that time. Images from that day were captured during measurement flight 20170614-1. Due to the sky conditions, flight altitude was limited to 300 ft (≈ 100 m), so no further altitude filtering was needed.

The measured ponds were assigned numbers 1 to 9 (Fig. 1). Pond nos. 1_2 and 4_2 merged with two larger ponds (nos. 1 and 4) over the course of 4 d. The selection of ponds for the analysis depended solely on the availability of in situ depth measurements carried out on-site as part of the measurement program and their location in the center of the flight pattern. Due to their position outside the photographed area, some ponds examined in König et al. (2020) had to be left out of this study.

Manual pond depth measurements were collected with a meterstick, and the locations of the measurements used here were determined relative to reference points beside the ponds or marked manually in aerial images from the previous day. We assume a horizontal accuracy of the measurement location of 0.3 m, which should not strongly impact the depth measurement accuracy in the center of ponds.

Pond depth measurements in the study area of PASCAL ranged from 7 to 26 cm on 10 June and 5 to 37.5 cm on 14 June (Table 1). We separated ponds by their either bright-blueish or dark-grayish appearance to prove the independence of our approach from the optical properties of the pond bottom. Manual depth measurements show that there were no systematic differences in depth between these groups. Drilling through the pond bottom after completing all other measurements revealed that all ponds were in an equilibrium state with sea level. Ponds 1, 3, and 4 were measured on both days, while ponds 5 to 9 were only measured on the first day. Ponds 1_2 and 4_2 had merged with no. 1 and no. 4 in the ongoing melt process. Pond no. 2 was measured only on the second day. Pond coverage increased notably in the study area during the 4 d.

Table 1 In situ pond depth measurement statistics and pond color type, either bright blueish (bb) or dark grayish (dg), from the PASCAL pond study site.

2.2 Aerial image data collected during MOSAiC 2020

Aerial image collection was part of helicopter grid surveys as part of the regular measurement program executed during the year-long drift campaign MOSAiC on board the RV Polarstern from 2019 to 2020 (Nicolaus et al., 2022). During the drift, RV Polarstern had to be repositioned several times, with a prolonged break in data during the initial pond formation period between 16 May and 19 June 2020 because of an inevitable crew exchange on Svalbard. At the time of RV Polarstern's return for expedition leg 4, a distinct floe had emerged that was round-shaped with a diameter of 1 km. It became the new location of the central observatory (in Webster et al., 2022, called Central Observatory 2, CO2). This floe of leg 4 formed from a former ice formation called the Fortress (von Albedyll et al., 2022), made of strongly compressed and deformed ice (second-year and multi-year ice, SYI and MYI) adjacent to the previous legs' Central Observatory area and surrounded by some FYI areas.

An orthomosaic and DEMs of the entire MOSAiC leg 4 floe are available from 30 June as well as 17 and 22 July 2020 (Neckel et al., 2022). They were compiled using the methods described in Neckel et al. (2023) and Fuchs (2023c) and cropped here to the extent of 2 km × 2 km around the floe. We classified the sea ice surface into three main surface types (ice/snow, ponds, and open water) by applying the sea ice image classification tool PASTA-ice (Proportional Analysis tool for Surface Types in Arctic sea Ice images) (Fuchs, 2023c) to the brightness-corrected orthomosaics (l2 data) on 17 and 22 July 2020 and to the brightness-corrected orthomosaic with cloud correction (l2b data) on 30 June 2020. The classification algorithm yields surface class maps in geospatial raster and vector data format, facilitating subsequent processing and analysis.

Photogrammetrically reconstructed DEMs from 30 June and 22 July 2020 had a vertical offset from zero caused in their processing. We leveled the open water level to zero using a flat plane fitted through all lateral snow/ice–open water boundary positions in the DEM within the cropped extent of 2 km × 2 km. These reference points were automatically extracted from the raster data DEM at the positions of touching surface class vector polygons. For 17 July 2020, this processing was not possible, as the DEM shows substantial deviations outside of the floe area because of the nonuniform drift of surrounding smaller floes during the grid survey, making it impossible to extract the water level from the DEM. However, manual inspection of level ice areas and single well-reconstructed water–ice boundaries confirmed that no further correction was required on this day in the MOSAiC floe area.

The floe edge was retraced manually in QGIS (QGIS Development Team, 2020) to retrieve statistical data only within the MOSAiC floe area.

3 Method development

The method development was performed with the PASCAL data due to the availability of ground-truth measurements. Flight patterns not yet adapted for depth determination also made it possible to develop a series of corrective measures, which, depending on the quality of the data, can also be used in future campaigns. All aerial images of a survey flight were taken with constant exposure settings and a mechanically and electrically fixed autofocus that was set to the flight altitude during pre-flight preparation.

3.1 Photogrammetric surface reconstruction

We use the commercial photogrammetry suite Agisoft Metashape to calibrate the camera optics and solve the complex aerial triangulation equations to calculate orthomosaics and DEMs as georeferenced raster data. The continuous drift of the ice during the measurement was thereby automatically corrected in the bundle-block adjustment by recalculating acquisition positions relative to the ice floe. Each surface point in the area of the studied ponds was captured on both days with more than nine different images. The ground sampling distance of both raster maps, orthomosaic and DEM, is 10 cm per pixel in the horizontal plane. The vertical resolution of the calculated DEM is 10 $\times {\mathrm{10}}^{-\mathrm{6}}$ m. Camera positions determined by aerial triangulation led to a reprojection error of images of 0.96 pixels (10 June 2017) and 1.03 pixels (14 June 2017). The ice drift was therefore successfully corrected by the determination of artificial image recording positions relative to the ice moved with the drift in a Lagrangian approach. Since ground control points (GCPs) with well-known position data were not available for a further accuracy assessment, a 2.00 m × 2.00 m black reference target located at the ice surface close to the ponds (Fig. 1) was used as a scaling reference. Length scale accuracy of the raster data was thereby determined to ± 2 %. DEMs were smoothed using the moderate depth map filter in Agisoft Metashape to avoid unnatural spikes in the model due to incorrect triangulation of single points.

Figure 2 Illustration of virtual and actual intersection lines in the photogrammetric reconstruction of melt ponds. Cameras 1 and 2 symbolize two individual pixels in images from different locations that captured the same point in a pond on the ice surface from a measurement angle α₁ and α₂. TP indicates the virtual point that is photogrammetrically reconstructed from the intersecting beamlines, assuming colinearity. AP indicates its actual position at the actual depth of the pond when refraction at the air–water interface is taken into account (change from angle of incidence β to the angle of emergence α). Δx and Δz quantify the horizontal and vertical deviation between the two, for which correction methods are introduced in the text.

3.2 Light refraction at the water–air interface

The photogrammetrical determination of the DEM relies on colinearity. Hence, optical beam paths between the observed ground surface and the lens are assumed to be straight. They are defined by the external orientation of the camera and the ground elevation. Distortions of the linear beam path are only considered in the camera optics. They are corrected with the Brown camera model integrated into the workflow. However, this basic assumption, valid for typical one-medium, low-level airborne observations, is invalidated in the case of underwater pond bottom observations by light refraction at the water–air interface. Due to the reduced speed of light in water compared to air, the electromagnetic wave and light beam are refracted more strongly away from the normal as they exit the water. The change from the angle of incidence β to the angle of emergence α at the interface between water and air is described by Snell's law:

$$\begin{array}{}\text{(1)}& n_{\mathrm{air}}\cdot \sin \left(\mathit{\alpha }\right)=n_{\mathrm{water}}\cdot \sin \left(\mathit{\beta }\right),\end{array}$$

with n representing the refractive indices of air and water. Both are assumed to be constant in this work with n_air=1 and n_water=1.335 (Millard and Seaver, 1990) (value for freshwater at 0 °C increased by 0.001 to account for salt remnants in the pond water). We found that dispersion, the wavelength dependency of the refractive index, does not have to be taken into account here in the wavelength range of the camera between 300 and 700 nm as it leads to deviations below the measurement resolution and accuracy. Further assumptions to describe the recorded beam paths from the pond bottom to the camera are the following.

1. There is no reflection or scattering of incident light at the pond surface or in the water column.

2. Reflection of incident light at the pond bottom can be approximated by Lambert's law.

These assumptions align with those made in Malinka et al. (2018) for the optical properties of pond bottoms with slight modifications. To match assumption 1, in clear-sky conditions, all images with sun glint on the pond surface need to be removed from the analysis. However, given the typically low solar elevations in the high latitudes, there was no sun glint on ponds in the observations used here due to (i) primarily wave-free pond surfaces and (ii) almost nadir measurements. The reflection of diffuse light at the pond surfaces in overcast conditions does not affect the structure recognition. It therefore does not pose a problem as it does in the pure optical retrieval algorithm of König et al. (2020).

Refraction at the pond–air interface may result in underestimating pond depth in photogrammetric measurements. In the following, we follow an idealized sketch of optical paths given in Fig. 2 to describe correction factors retrieved to cope with the impact of refraction. The overarching goal of the method was to preserve or restore colinearity in the multi-view surface reconstruction so that an integrated evaluation with Agisoft Metashape is still feasible. This approach therefore differs strongly from previous studies, which were set up on a completely new and complex set of equations (e.g., Westaway et al., 2001) and then take into account the refraction of light but do not rely on highly specialized programs to solve it efficiently or merely neglected the impact of optical effects (e.g., Casella et al., 2017; Divine et al., 2016).

Figure 3 Geometrical evaluation of Eq. 2 (including Eq. A6), showing how the specific horizontal mismatch κ changes with two measurement angles α₁ and α₂. We found that below 40° measurement angle errors induced by the horizontal mismatch remain small enough to be neglected in the automatized reconstruction.

Figure 2 shows the virtual and actual intersection of beamlines in a pond from two opposite and monocular observational positions. The actual intersection point AP represents the true pond bottom, while the virtual point TP is located at the depth where the beamlines would intersect without refraction. The latter is the one determined by Agisoft Metashape. To better understand the approximations derived from it, we retrace the optical path from the pond bottom to the camera in three steps: (I) a recognizable pattern allows a point on the pond bottom to be clearly identified on different images captured from different positions. We call this key point at its actual position AP and assign the Cartesian coordinates $\mathrm{AP}(X_{\mathrm{AP}},Y_{\mathrm{AP}},Z_{\mathrm{AP}})$. (II) As we consider AP a Lambertian reflector, identical beams radiate from the point in all unobscured directions. Two of them reach the monocular observation points camera 1 and camera 2, which symbolize two pixels in different aerial images taken along the flight track. These beamlines are not straight but bent at the water–ice interface. The angle of emergence α is larger than the angle of incidence β and defined by Snell's law: $\mathit{\alpha }=\arcsin (n_{\mathrm{water}}\cdot \sin (\mathit{\beta }\left)\right)$ (Eq. 1 rearranged). (III) The coordinates of the origin of AP are determined by aerial triangulation from all available camera positions. Since this is based on the assumption of colinearity, the obtained virtual position of the point in the reconstruction is located in $\mathrm{TP}(X_{\mathrm{TP}},Y_{\mathrm{TP}},Z_{\mathrm{TP}})$. TP deviates in height from the original position Z_AP by Δz and since camera positions usually do not all have precisely the same elevation angle also in its horizontal position X_TP and Y_TP by Δx. This results in two deviations in the colinearity approximation caused by the pond water that must be corrected or avoided, a vertical and a horizontal one. Both deviations are discussed separately in the following two paragraphs.

Figure 4 Geometrical evaluation of depth correction factor γ for equal measuring angles (a) and differing measuring angles (b). n_w is the refractive index of water. The vertical axis of panel (a) is split at γ=1.527, which equals the depth correction factor for the limited measurement angle α=40° (indicated by the black lines). γ strongly increases for α>40°. In panel (b), the 40° threshold is marked with the purple lines.

3.2.1 Horizontal deviation

The horizontal deviation Δx potentially causes a mismatch of point detections and should therefore be avoided in an integrated scheme. We consider it sufficiently suppressed when Δx is smaller than the ground sampling distance. Δx is directly dependent on both angles of emergence and the measured pond depth Z_TP. We define the ratio between horizontal deviation and measured pond depth as deviation factor specific horizontal mismatch κ, with

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{\kappa }={\displaystyle \frac{\mathrm{\Delta }x}{Z_{\mathrm{TP}}}}.\end{array}$$

Figure 3 shows how κ changes with both incident angles. At a flight altitude of 300 ft, which is a good choice for highly resolved pond studies, the ground sampling distance of the CANON camera system is approximately 0.05 m. This means that for all measured pond depths Z_TP up to 1.5 m, the maximal horizontal deviation $\mathrm{\Delta }x=\mathit{\kappa }\cdot Z_{\mathrm{TP}}=$0.042 m caused by refraction remains below the measurement resolution and the photogrammetric projection error when we restrict incident angles to α_max=40°. We see it as a good compromise between moderate horizontal deviation and enough field of view to preserve sufficient overlap of images. Therefore, all image pixels at larger opening angles relative to the nadir are neglected in the reconstruction. This is done by creating masks individually for every single image depending on the orientation angle of the camera derived from the camera alignment process in Metashape.

Figure 5 Illustration of different elevation and depth levels relevant for the determination of individual pond depths. Schematically, h_surf is the vertical position of a pond surface, h_bot is the vertical position of its pond bottom, and h_sea is the vertical position of the sea surface. The bathymetric map B(X,Y) describes the actual pond depth at the projected geographic north and east coordinates (Y and X) in ponds that are not uniformly deep in reality. T(X,Y) is the height of the pond margin at coordinates X and Y.

3.2.2 Vertical correction

After examining how horizontal deviations can be avoided, this section concerns the underestimation of the measured pond depth Z_TP owing to refraction. We discuss how it can be corrected with a correction factor γ defined by $Z_{\mathrm{AP}}=\mathit{\gamma }\cdot Z_{\mathrm{TP}}$. γ(α) is given by

$$\begin{array}{}\text{(3)}& \mathit{\gamma }\left(\mathit{\alpha }\right)=\cos \left(\mathit{\alpha }\right)\cdot \sqrt{n_{\mathrm{water}}^{\mathrm{2}}-{\sin }^{\mathrm{2}}\left(\mathit{\alpha }\right)}.\end{array}$$

It can be shown that the correction factor γ converges towards the refractive index of water n_water for arbitrarily small α and eventually becomes equal to n_water at exactly α=0 (see Appendix A for the equation solutions).

$$\begin{array}{}\text{(4)}& \underset{\mathit{\alpha }\to \mathrm{0}}{lim}\mathit{\gamma }\left(\mathit{\alpha }\right)=n_{\mathrm{water}}\end{array}$$

This contradicts the statement in Casella et al. (2017) that refraction does not influence underwater depths retrieved from almost nadir images. Instead, our mathematical and geometrical evaluation shows that almost nadir depth measurements must be multiplied by the refraction index to achieve correct depth values. When α increases, γ rises slowly at first and eventually very strongly (Fig. 4a). In the previously defined range of maximal 40° incident angle, γ reaches a maximum value of γ_max=1.527. However, since derived depths result from averaging numerous intersection rays with mostly small angles of incidence, the small increase in γ is ignored, and γ is kept constant and equal to the refractive index of water n_w in our method for the correction of all underwater pixels.

Wright and Polashenski (2018)Miao et al. (2015)Miao et al. (2015)Miao et al. (2015)

Table 2 Pixel-wise input features to the classification scheme.

An additional conclusion can be drawn for the horizontal deviation from this analysis. Figure 4 shows that water depths are systematically underestimated with the measured pond depth, which means that the previously obtained maximum opening angle α_max in fact allows for greater actual pond depths than assumed in Sect. 3.2.1, in which we limited the virtual pond depth Z_TP to 1.5 m.

Table 3 Sea ice surface type classes used in the classification tool PASTA-ice to semantically segment orthomosaics.

3.3 Pond depth determination

Pond depth d is the vertical extent of the water column in ponds. It is composed of the vertical position of the pond bottom h_bot and the height of the pond water surface h_surf (Fig. 5). Large altitude inaccuracies of aircraft positioning systems, especially at high latitudes, make it impossible to use the GPS aircraft altitude and the reconstructed distance to the ground as an absolute height reference above sea level. That is why studies of ice topography with, for example, ALS usually use areas of very thin ice or open water as a reference to water level (e.g., Hutter et al., 2023). For ponds, the reference is even more complex since pond depth is, as previously mentioned, not only prescribed by the topography of the pond bottom but also by the individual height of the water level in the pond. This water level is typically above sea level, partially caused by impermeable sea ice or later in the season, when ice is typically permeable, by the density difference between freshwater in the pond and underlying seawater. The pond water level is individual for every single pond, especially during the early stages of melt pond formation, i.e., the time in melt season before the ice gets permeable enough to allow for pond drainage (see stage I, Eicken et al., 2002). Hence, we calculated relative pond depth without an absolute reference. To this end, the water level in each pond was determined by its height in the DEM at the edge of the pond. The method therefore highly benefits from optical aerial images enabling high-resolution surface type classification with PASTA-ice (described in the next section) on exactly the same raster as the reconstructed topography data. We derived the pond margins as vectors (line strings) from the classified images and overlaid the photogrammetrically retrieved DEM raster data T(X,Y) with the pond margins. Then, we extracted the relative height of the pond surface h_surf at the pond margin from the DEM, where this marks the transition from ice to pond in the smoothed topography. Due to the smoothing, we expect the method to also work for ponds with almost vertical walls later in the season (Fetterer and Untersteiner, 1998), which were not part of the evaluation set. This extraction was done using the Python libraries geopandas, rasterio, and rasterstats (Jordahl et al., 2020; Perry, 2015; Gillies, 2013, respectively). The two-dimensional bathymetric map B(X,Y)_i of each individual pond i, with X and Y as geographic north and east coordinates in the projected coordinate system, was subsequently retrieved from

$$\begin{array}{}\text{(5)}& B(X,Y{)}_i=-(h_{{\mathrm{surf}}_i}-T(X,Y))\cdot n_{\mathrm{water}},\end{array}$$

where n_water is the refractive index of water as a depth correction factor as discussed in Sect. 3.2.2. Depth is specified positive down. Last, we compiled a pond-depth-corrected topography map T_pond(X,Y) from the individual pond bathymetries.

3.4 Pond classification with PASTA-ice

For the automatic detection of ponds in aerial images, we used the Proportional Analysis tool for Surface Types in Arctic sea Ice images (PASTA-ice, source code repository in Fuchs, 2023a, previously used in Thielke et al., 2023; Niehaus et al., 2023). In the following, we briefly summarize the algorithm due to its importance for workflow automation. Details, including the method development and evaluation, are found in Fuchs (2023c).

PASTA-ice is tailored to the aerial images captured with the AWI imaging system. Besides focusing on the semantic separation into surface type classes, the algorithm aimed at retracing pond outlines used to extract pond levels and to compile statistics on individual ponds. Image classification is done pixel-wise in brightness-corrected orthomosaics (e.g., Neckel et al., 2023) based on absolute R, G, and B values and ratios thereof (Table 2).

Classification is performed with the random forest classifier implementation in Scikit-learn (Pedregosa et al., 2011). The classifier was trained and tested using data from manually selected areas in very diverse sea ice surface appearances recorded during PASCAL. Pixels in orthomosaics are classified into nine different sea ice surface subclasses that belong to three main classes (Table 3): snow/ice, open water, and ponds (including submerged ice). Adjacent subclass pixels of similar main classes are subsequently combined into main class vector objects if these consist of, at minimum, 100 pixels (the threshold was chosen similar to Huang et al., 2016). This threshold is applied as a minimum area requirement to match the baseline of high-resolution data that objects are resolved from various pixels (e.g., Wright and Polashenski, 2018). The chosen threshold corresponds in the orthomosaics of this study to an area of approximately 1 m² (PASCAL) and 25 m² (MOSAiC). Smaller objects are considered noise and are added to the largest adjacent object; their area fraction is taken into account when estimating the inaccuracy of the classification result.

For spatial analysis, main class objects are converted to polygon geometries defining the outer and, if present, inner edges of the object. They also include an attribute table that contains information on the subclass proportions and a classification confidence proxy from the prediction probability output of the classifier (Fuchs, 2023c). Classification recall and precision are somewhat limited for the very specific sea ice surface subclasses listed in Table 3 due to overlaps in appearance but high for the combined main classes (Fig. 6). High accuracy values mainly result from large sample sizes, resulting in large numbers of true negative pixels. All pond objects are reclassified to submerged ice if they are located spatially between a snow/ice and open water object and if the open water object is larger than the pond object.

The polygons of ponds, i.e., their outlines defined by vertices and connecting edges that resemble the snow/ice to pond interface, are used to extract the pond level from the DEM. The very high accuracy in classifying the main classes indicates a sufficient detection of the transition areas from pond to ice and, thus, of the pond polygons.

Figure 6 Accuracy, recall, and precision evaluation for subclasses and main surface type classes in the PASTA-ice classification scheme retrieved from a test dataset compiled from images collected during PASCAL. The orange and green lines mark the separation between open water and ice classes and between ice and pond classes, respectively. Detail of Fig. 3.19 (c, f) in Fuchs (2023c).

3.5 Pond margin detection and height correction in the PASCAL data

Pond margins in the PASCAL study area were traced manually in orthomosaics using QGIS (QGIS Development Team, 2020) to better assess the depth correction algorithm without the impact of any classification inaccuracies expected in this particular study area (traced polygons shown in the study site overview, Fig. 1). Due to the deformed ice surrounding this specific location, shadows tended to impact the automatic classification scheme in PASTA-ice on this small scale, which could eventually strongly falsify the pond exterior detection needed to derive h_surf, especially since misclassified shadows that stretch from ponds into adjacent ridges can be partially far above the water level. In larger sample sizes and more even ice areas, where most ponds usually form during summer melt, automatic surface classification with pond margin detection is easily possible, as shown in Sect. 5 with the MOSAiC data.

Figure 7 Photogrammetrically reconstructed DEM of the study site on 14 June 2017 with marked ponds (no. 1–no. 4 from Fig. 1). Missing GCPs led to a strong tilt of the surface. Inaccurate GPS heights and reconstruction in the WGS84 ellipsoid cause an absolute offset from 0 m.

Missing GCPs and slightly arbitrary camera optics mainly caused by a bubble-shaped protection window in front of the moveable mounted camera – for shock protection – furthermore led to large-scale deviations like curved and tilted surfaces in the Agisoft DEM retrieval of the PASCAL data (Fig. 7). These deviations were approximated as a linear slope on the length scales of ponds. As a correction, we fitted a two-dimensional plane through each pond outline in the DEM. This regressed plane is assumed to represent the water level $h_{{\mathrm{surf}}_i}$ in Eq. (5).

$$\begin{array}{}\text{(6)}& h_{{\mathrm{surf}}_i}(X,Y)=a\cdot X+b\cdot Y+Z\end{array}$$

Z is the absolute height correction resulting from the mean vertical deviation of the retrieved surface from sea level. This mean vertical deviation is mainly caused by inaccurate GPS altitudes at higher latitudes. Height levels along the corrected pond edge deviate only slightly from zero (mostly less than ± 0.05 m) as shown in the box-and-whisker plots in Fig. 8. In bending-free DEMs (e.g., the ones from MOSAiC presented below), the reference height can be derived solely from the mean elevation of the pond outline $h_{{\mathrm{surf}}_i}$.

Figure 8 Deviation of the pond outline topography given by T(X,Y) to the fitted pond surface planes $h_{{\mathrm{surf}}_i}$ on 10 June 2017 (a) and 14 June 2017 (b).

4 Results of the method development

4.1 Photogrammetrically derived pond depth from PASCAL compared to manual measurements

Bathymetric charts of ponds in the study area were calculated from aerial images for all manually sampled ponds as previously described (Fig. 9). To account for the location inaccuracy of the in situ data, photogrammetrically derived pond depth data were averaged in a circle with a radius of 0.3 m around the point measurements (Fig.  9).

Figure 9 Photogrammetrically derived pond bathymetry charts of all ponds in the PASCAL study area, for which in situ data (dots) were available. The colored shapes show the photogrammetrically reconstructed pond bathymetries and the dots show the in situ measurements. The size of the dots equals the area used for averaging the photogrammetrically derived depths before comparing them to the in situ data.

In situ and photogrammetry yielded almost identical pond depths (Fig. 10). The original pond color (see Fig. 1) indicated by the color of the dots (blueish, bright, or grayish dark ponds) in Fig. 10 does not affect the reconstruction. However, pond size shown by the size of the dots in Fig. 10 seems to have one specific influence on the reconstruction: small ponds (<1 m in diameter), indicated by tiny dots, are underestimated in their pond depth. As a bias of the measurement method, we retrieve bias $=-\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{3}}$ m with a root mean square error of RMSE $=\mathrm{3.84}\times {\mathrm{10}}^{-\mathrm{2}}$ m and mean absolute error of MAE $=\mathrm{2.65}\times {\mathrm{10}}^{-\mathrm{2}}$ m.

Figure 10 Pond depth comparison between manual in situ data and photogrammetrically derived bathymetry. The dot colors indicate the original pond color (blueish, bright, or grayish dark ponds); see Fig. 1. The size of the dots indicates melt pond size ranging from 0.5 to 116.9 m².

4.2 Photogrammetrically derived pond depth from MOSAiC compared to echo-sounding data

It was noted above that ridge structures directly adjacent to the ponds in the PASCAL study area required a manual tracing of pond polygons to get a precise reference water level. Improved flight patterns during the MOSAiC expedition (mowing-the-lawn pattern) with a regular lateral and forward overlap in images, together with well-classifiable surface conditions, increased the capability of the algorithm to work entirely autonomously. We evaluate the photogrammetrically derived pond depths with the newly developed autonomous pond investigation system Böötle, equipped with a downward-facing echo sounder (Oppelt and Linhardt, 2023). We compare pond depths of a lake-like pond, Mystery Lake, which reached depths of more than 2.5 m. It was regularly mapped during helicopter survey flights and with Böötle. We chose the two datasets collected closest in time on 7 July 2020 (helicopter survey) and 9 July 2020 (echo-sounder measurements). However, the temporal difference of 2 d restricts us from retrieving precise errors and corrections, but we can still use it to confirm the overall method. The comparison is particularly interesting because the pond exceeded the maximum depths assumed in the method development, and due to smaller lateral image overlap, the opening angle in the acquired images could not be limited to 40°. Figure 11 compares the photogrammetrically derived and echo-sounding pond depths. Overall, both methods agree in pond depths for Mystery Lake. Yet, the data show a divergence at greater depths (>1 m) that can be attributed to either the further deepening of the pond within the 2 d, a systematic underestimation of the reconstructed depth, or both. The latter would be a consequence of the unrestricted opening angle, which would actually require a larger correction factor than what is applied. The percentage deviation was −11.5 % at depths below 1 m and increased to −24.2 % in greater depths.

4.3 Impact of flight pattern

Since the mowing-the-lawn pattern is among the most time-consuming flight patterns, we further determined to what extent other flight patterns, for example, straight flight legs with sufficient forward overlap, also lead to reasonable results. Based on pond no. 1 (10 June 2017), we investigated how the measurement accuracy depended on the overlap of images and the lateral offset of the flight lines. Figure 12 shows bathymetric charts of pond no. 1 retrieved from (a) a few measurement positions along a straight line (4 images), (b) more measurement positions along a straight line (10 images), (c) similar measurement positions with lateral offset (9 images), and (d) many measurement positions with lateral offset and all available measurement positions as used for the comparison with in situ observations before (31 images). The camera matrix and image recording positions were optimized before and kept constant during this study so as not to impact the error.

Figure 11 Hexagonal binning diagram showing the comparison of pond depth measurements in Mystery Lake during MOSAiC from an in situ echo sounder (deployed on the platform Böötle, 9 July 2020) and photogrammetrically derived bathymetry (helicopter survey flight, 7 July 2020). A total of 7186 measurement points were compared. Only bins that contain at least one point are shown. Percentage deviations (PDevs) are calculated for echo-sounder depths below and above 1 m (separated by the dashed red line).

Figure 12 Photogrammetric pond bathymetries of pond no. 1 (10 June 2017) derived from different subsets of image recording positions used in the reconstruction. (a) Reconstructed from 4 images along a straight line with almost nadir measurements. (b) Reconstructed from 10 images along a straight line. (c) Reconstructed from 9 images scattered above the pond with almost nadir measurements. (d) Reconstructed from 31 images scattered all around the pond. RMSEs are calculated by comparing the reconstructions to the in situ measurements.

The accuracy strongly increases with an increasing number of measurement positions, larger incident angles, and especially lateral offset of the measurement positions when comparing measurement positions on one straight line or isotropically distributed over the pond (Fig. 13). A reconstruction from 10 images recorded along a straight line yields a higher error of 7.90 cm compared to 9 images with both lateral and forward overlap (6.98 cm). The quantitative differences are still small, while subjectively, the differences in Fig. 13 even exceed these error estimates. Lateral offset can best be achieved with a mowing-the-lawn flight pattern.

5 Melt ponds on the MOSAiC floe

Regularly flown floe grids during MOSAiC (with mowing-the-lawn flight pattern) combined with the newly developed pond bathymetry retrieval enable an unprecedented three-dimensional analysis of melt ponds. With the flight pattern, lateral and forward overlap in images was achieved, and owing to open water areas around the major floe and the DEM correction with ALS data by Neckel et al. (2023), a leveling of all data to sea level was possible. To do so, we fitted a flat plane through all snow/ice–open water edge positions in the DEM around the floe and subtracted it from the DEM. For the first time, we could thus retrieve pond bathymetries, relate pond level to sea surface height, and track pond changes from flight to flight with aerial imaging.

5.1 Evolution of pond coverage, bathymetry, level, and volume

The first pond formation on the MOSAiC floe happened in May 2020 as observed from satellite images (Webster et al., 2022), at the time of the absence of RV Polarstern from the MOSAiC floe for crew exchange. After the return in mid-June, continuous pond formation was observed. Pond coverage on the MOSAiC floe between the end of June and late July varied between 22.0 % and 23.7 % (Table 4).

At the end of June, the pond cover was well developed and amounted to 22.3 % (Fig. 14). By mid-July (11–13), most ponds had drained (Webster et al., 2022). Having orthomosaics of the entire floe from 30 June 2020, 17 July 2020, and 22 July 2020 (Fig. 14), we thus lack comprehensive aerial data on the days of drainage but cover the period both before and after. On 30 June, several vast, exceptionally deep ponds (>2 m) had formed on the MOSAiC floe, along with many smaller ponds (Fig. 14). After drainage, the pond cover became more fragmented and braided-like, with more shallower ponds (<1 m) in the center of the floe. We discuss this significant loss in depth of the largest ponds (including Mystery Lake) in the Discussion section. A direct comparison of the conditions before and after the drainage event shows that individual ponds in the strongly deformed center of the floe and smaller ponds have remained unchanged in their shape or have even grown. All prominent large ponds on the floe underwent major changes during drainage and show bare patches of ice that were formerly submerged under pond water. Nevertheless, total pond coverage was relatively constant at 22 % at the time of the three measurement flights (Table 4). A possible underestimation of the pond coverage from the PASTA-ice classification was relatively strong on 30 June and 22 July 2020 at 3.7 % and 3.9 %, respectively. This possible underestimation was presumably caused by small pond objects in the classification output that were below the 100-pixel minimum threshold for objects to be resolved, as described in the “Method development” section.

Previous studies found that pond surfaces before drainage were above sea level, forming a hydrostatic head (e.g., Perovich et al., 2021, based on SHEBA data). On the MOSAiC floe, on 30 June, 2 weeks before the main vertical drainage event, we find, in contrast, that 90.1 % of the total pond area was already close to sea level (less than 0.2 m above sea level). Only a few small ponds (7.1 % of the entire pond area) in the strongly deformed center of the floe were well above sea level (more than 0.3 m above sea level). The 50 largest ponds (covering 64.7 % of the total pond area) were closer to sea level and therefore had a low hydrostatic head.

The partitioning of meltwater has raised particular interest in MOSAiC since distinct freshwater lenses were observed around and below the floe, strongly impacting the physical and biological sea ice system (Smith et al., 2023). Melt ponds form a reservoir in the meltwater budget. Their bulk volume is therefore of great interest to better assess and understand the total budget. Having derived pond bathymetry maps of the entire floe, we can, for the first time, derive the overall volume of this meltwater reservoir directly, in contrast to extrapolating it from single transect lines, which was the only available method before (e.g., Perovich et al., 2021; Webster et al., 2022). The bulk volume of meltwater in ponds on the floe results from the pond bathymetry and the area covered by ponds. Deriving it as an area-specific value, we find that area-specific pond volume (asPV) on the entire floe was fairly constant at 0.040 to 0.051 m³ m⁻² (volume meltwater per floe area, Fig. 15a).

The floe-wide, high-resolution data give us the opportunity to further investigate the spatial and temporal evolution of the pond volume. Before drainage, half of the meltwater in ponds was stored in particularly wide ponds with large diameters, quantified here by the diameter of the biggest disk that can be fitted into the pond shape. With the drainage event, such wide ponds like Mystery Lake disappeared (Fig. B1e), and ponds with diameters 5 to 20 m got the largest sinks for meltwater (92.5 %). Some of them were particularly large in area (up to 13 000 m²) but simultaneously small in diameter (<15 m) (Fig. 15b, pond marked in dark khaki). After the drainage event, only half of the meltwater (17 July 2020: 53.6 %, 22 July 2020: 52.7 %) was stored in ponds with diameter >12.5 m. Ponds of similar diameter increased in depth with time by 2.5 cm between each flight (median depth increase) and therefore compensated for the loss of volume in the particularly wide ponds that vanished with the drainage event and, with that, decreased the mean pond depth (mD) by up to 4 cm (Fig. 15c). At the same time, the overall melt pond coverage (MPF) remained relatively constant. Still, ponds changed their appearance strongly to braided-like pond patterns with smaller diameters but large connected pond areas (Fig. 14). It is generally noticeable that the pond depth correlates strongly with the diameter size introduced here, with larger diameters allowing for greater depths (Fig. 15c). Following the shape of the depth d to diameter ϕ dependency, we approximated a square root function with the least square method to the binned data (min. three ponds per bin). For the time before the drainage event, we yield $d\left[m\right]=\mathrm{5.36}\times {\mathrm{10}}^{-\mathrm{2}}\sqrt{\mathit{\varphi }\left[m\right]}$. After drainage, the relation changes to $d\left[m\right]=\mathrm{6.61}\times {\mathrm{10}}^{-\mathrm{2}}\sqrt{\mathit{\varphi }\left[m\right]}$. Both fits approximate the individual pond depths with an RMSE of 5.5 to 5.7 cm. However, the continuous deepening of the ponds after the drainage event clearly shows that this fit is limited in time. The shallower depth of smaller ponds explains why they contribute less to the total pond volume than to the total surface area of ponds, visible through the slightly skewed distribution function (Fig. 15b). We also included pond volume fraction, areal fraction, and depth, resolved over pond area in Appendix B (Fig. B1). However, signals relative to the pond area are much less pronounced than if they are separated by the diameter measure used here.

Figure 13 Comparison of the reconstruction accuracy using images only along a line or from horizontally distributed positions. RMSE of pond no. 1 depths (10 June 2017) derived from the comparison of in situ measurement points to reconstructed pond depths with the same logic as in Fig. 10. Shown are RMSEs for different subsets of image recording positions used as input to the reconstruction. Differentiation is made between points along a single line and points scattered isotropically above the pond. The smaller the smallest measurement angle, the farther away the most distant recording point is and the larger the offset.

Table 4 Pond coverage on the MOSAiC leg 4 floe retrieved from helicopter aerial imaging and PASTA-ice classification on 3 different days. Numbers in brackets depict the confidence range obtained from a newly developed confidence evaluation (Sect. 3.4).

Figure 14 Overlays of photogrammetrically reconstructed pond data on orthomosaics of the MOSAiC floe from different MOSAiC floe grid surveys (orthomosaics and DEM from Neckel et al., 2022, reprojected in UTM31N). The panels show pond bathymetry (a) and pond level above sea surface height (b) reconstructed from an airborne survey flight on 30 June 2020, as well as pond bathymetries from flights on 17 July 2020 (c) and 22 July 2020 (d). Ponds were detected automatically using the PASTA-ice surface type classification.

Figure 15 Distribution functions of pond volume (a), pond areal fraction (c), and mean pond depth (b) in 13 pond diameter size bins from 0 to 30 m and larger than 30 m on the MOSAiC floe. Data retrieved from photogrammetric reconstructions from aerial images collected on survey flights on 30 June (pre-drainage), 17 July, and 22 July 2020 (post-drainage) are shown. The total floe-area-specific pond volume (asPV) is given in the legend of (a), and the total pond coverage (MPF) on the floe is in the legend of (c). Mean pond depths (mD) on the entire floe are listed in (b). Vertical error bars in (b) depict the standard deviation of pond depth in the diameter bin. For better visibility, their x position in the graph is slightly offset. Fits show square root functions fitted to the binned data with the least square method. Filled areas in (a) and (c) show color-coded contributions of the two ponds, whose bathymetry and shape are shown on the lower right. The color map of the bathymetry reaches from 0 to 2 m depth (Fig. 14). Red circles show the largest disk that can be fitted into the pond shape, which is used here to quantify pond diameters.

5.2 Comparison to other in situ and satellite observations

Measurements on MOSAiC were carried out using many different methods and may lead to different results. In the following, we compare the aerial derived data to available results from high-resolution satellite observations (Webster et al., 2022, using the Wright and Polashenski, 2018, classification algorithm OSSP) and in situ transect lines (Webster et al., 2022) to assess the accuracy of our results and the representativeness of observed areas. Transect lines were repeatedly revisited paths on which extensive in situ measurements of ice thickness, snow, and ponds were carried out. The transect considered here was the longest and surrounded the entire floe (Webster et al., 2022, Fig. 2). To compare to these transects, we derived pond properties additionally within a 10 m buffer zone around retraced transect footpaths in the aerial images. Satellite and aerial-derived pond coverage of the floe is temporally sparse, but both remain relatively stable at around 22 % over the entire observational period with differences <2 %. Along the transect lines, helicopter-derived pond coverage increases from ∼ 10 % in late June to ∼ 20 % in mid-July. It thus resembles the evolution of the in situ observed coverage (Fig. 16a). Small under- and overestimations of 2 % to 5 % occur, probably caused by collocation inaccuracies. All products thus seem to be comparable for their observed area. However, it becomes apparent that along the floe edge, pond coverage doubled between the end of June and mid-July and thus acted differently from the relatively constant pond coverage on the entire floe.

Pond depth along the transect line was observed to be about 2.5 cm higher (∼ 15 % to 20 %) than the mean derived pond depth for the same area from the aerial-derived bathymetric maps (Fig. 16b). The (very variable) pond depth on the entire floe is underestimated by the transect area in June. It matches well in July after drainage, probably because the extraordinarily deep ponds in the floe center flattened, making the pond cover more uniform. Both methods resolved a slight increase in pond depth at the floe edge. Before harmonization through the drainage event, the very deep ponds in the middle of the floe, a region not covered by the transect studies, caused a more constant pond volume than extrapolated by Webster et al. (2022) from the transect area (Fig. 16c). Good agreement between the two methods along the transect lines suggests that differences occur mainly due to different observed areas instead of methodical differences.

Figure 16 Time series of melt pond coverage (a), mean pond depth (b), and specific pond volume (c) retrieved from in situ transects, high-resolution satellite images (both Webster et al., 2022), and helicopter-borne aerial images. Floe properties are retrieved for (aerial and satellite) or extrapolated to (in situ) the entire MOSAiC leg 4 floe. Transect properties are measured along the transect line (in situ) or retrieved within a 10 m buffer zone around the retraced line (aerial). The in situ specific pond volume is the pond volume presented in Webster et al. (2022) divided by the floe size. However, due to a correction, we recalculated the specific pond volume using the mean pond depth and pond coverage values (in consultation with the authors). Both lines are included to avoid misunderstandings (see Webster et al., 2024).

5.3 Upscaling factor for in situ depth measurements

In situ measurements in ponds are often restricted to very few single points. We investigate whether and how the mean pond depth of entire ponds can be extrapolated from single in situ point measurements by subsampling our high-resolution photogrammetric pond bathymetry reconstruction. To this end, we benefit from the unprecedented dataset in resolution within ponds and the total number of ponds.

We subsampled pond depths from survey flights on 30 June (pre-drainage) and 22 July 2020 (post-drainage). For each pond, we extracted the point furthest away from the pond edge as the pond center, the so-called pole of inaccessibility (PIA). We assumed that in situ depth measurements are typically or can be easily performed at this point, expecting the most representative depth in the center. All pond objects on the major floe classified with PASTA-ice were considered. Except to prevent errors caused by smoothing of the DEM in small ponds, ponds with a maximum distance between the center and pond edge of <1 m were neglected. In total, the evaluation was based on 1.6 ×10⁶ pixels in 1621 aerially observed ponds. A selection between older (second- and multi-year ice) and younger ice (first-year ice) was based on personal testimonies reporting a younger ice area in the south of the floe (orientation of the floe in June–July 2020) and older ice in the center.

Comparing mean pond depth d_mean derived from the entire pond bathymetries with single pond depth measurements in the center d_center reveals that a single measurement in the center of the pond strongly overestimates the mean depth. On average, d_mean is only 52 % of d_center (Fig. 17). Hence, the mean pond depth and pond volume are much smaller than assumed based on single measurements from the pond's center. A descriptive form factor κ, which we define as

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathit{\kappa }{\displaystyle }={\displaystyle \frac{d_{\mathrm{mean}}}{d_{\mathrm{center}}}}=\mathrm{0.52}\end{array}$$

and retrieve from

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathit{\kappa }{\displaystyle }={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\left({\displaystyle \frac{d_{\mathrm{mean}}}{d_{\mathrm{center}}}}\right)}_i,\end{array}$$

including all n=1621 ponds, shows strong statistical significance in the large dataset. Residuals of the linear fit and thus deviations from this generalization are normally distributed (Fig. 18), independently if all ponds are considered (Fig. 18b) or if subsets are taken for ponds in the pre- and post-drainage state (Fig. 18f) on younger less deformed (FYI) or older strongly deformed (MYI) ice (Fig. 18d). Furthermore, residuals do not correlate with the pond area (Fig. 18a), the maximum disk area from which the single depth measurements are taken (Fig. 18c), the pound roundness (defined here as $\frac{\text{max. disk area}}{\text{pond area}}$) (Fig. 18e), the single depth measurement (Fig. 18b), or the drainage state (Fig. 18f). Only ponds on young ice seem to develop a slight dependency of the form factor fit residuals from the pond depth with a Pearson correlation coefficient of −0.41 (Fig. 18d). From this robust relationship between d_center and d_mean independent of the pond type, we conclude that the descriptive form factor κ=0.52 can be used to upscale single point measurements from the center of the pond (pole of inaccessibility) to mean pond depth.

The introduced form factor describes the average relationship between center depth and mean depth. For individual ponds, deviations may occur as indicated by the residuals of the fit. We explore the number of sampled ponds at which the form factor becomes a valid estimator for the mean pond depth of multiple ponds. We calculated the form factor for 100 randomly chosen subsamples of different sizes n from the 1621 ponds, similar to bootstrapping. The defined goal was that 99.3 % of the retrieved form factors within the 100 random samples are located within the range of  ± 10 % of the form factor κ=0.52. Figure 19 shows the box-and-whisker plots for different pond sample sizes n, retrieved from the 100 random samples per size n. The upper and lower whiskers define the 99.3 % range. It can be seen that from 64 ponds on, the requirement is fulfilled. Thus, at least 64 ponds are needed to yield a valid approximation of the mean pond depths from single point measurements in the center of the ponds, applying the form factor κ.

6 Discussion

6.1 Photogrammetric reconstruction workflow for melt pond bathymetries

Photogrammetry can be used to reconstruct the sea ice surface topography from aerial images (Neckel et al., 2023; Divine et al., 2016). In the first part of this study, we investigated the ability to also reconstruct the full melt pond bathymetry. In particular, we have developed methods to consider water in the light paths used for reconstruction, in addition to air. We developed the methods using data from the PASCAL campaign. For this particular campaign and study site, more extensive processing was necessary, including corrections due to curved DEMs and the manual digitization of the pond edges. The evaluation of MOSAiC data showed that these were no longer required with better data quality through adapted flight patterns and better-classifiable images. We thus present an algorithm that can be adjusted depending on the quality and method of the recorded data. In the following, we discuss the results obtained for the bathymetry reconstruction.

6.1.1 Impact of sky conditions, pond color, size, and shape

Within the presented dataset, no difference in the accuracy of the depth retrieval could be traced back to sky conditions or pond colors. However, we identified potentially larger uncertainties for very small and large ponds. Ponds smaller than 1 m in the horizontal expansion (equalling <10 times the ground sampling distance) were not properly reconstructible. We hypothesize that the essential smoothing of the image depth maps flattens the topography in such small ponds. Further, on very large ponds, we assume that wind fetch causing ripples on the water surface could lead to disturbing reflections. In a weak form, this effect was observed on the abovementioned Mystery Lake (diameter approx. 80 m) on the MOSAiC floe, but it was not strong enough to degrade the retrieval quality.

In addition, unavoidable smoothing of the topography in the photogrammetric calculations may also lead to larger uncertainties for steep pond bottom slopes at the outer edges of ponds. Too few measurements were available to evaluate this quantitatively. However, we observed such effects by comparing the results visually to data presented in König et al. (2020). Since their method relies on spectral differences, it does not smooth the bathymetry toward the pond edges. Due to a lack of data, we could not test whether vertical pond walls that occur later in the season negatively impact the pond surface extraction from the DEM. However, we assume that, in this case, the automatic smoothing reduces the error.

Figure 17 Linear dependency between single pond depth measurements at the pond center and mean pond depth reconstructed from aerial images of the MOSAiC leg 4 floe on 30 June 2020 (pre-drainage) and 22 July 2020 (post-drainage). Each dot represents a single pond. Particularly deep ponds were circled to make the dots more visible. All data were retrieved from aerial observations.

Figure 18 Residual analysis of the linear fit given in Fig. 17. Residuals of all ponds are shown with respect to their area (a), center depth (b), pole of inaccessibility (PIA) disk area (c), center depth and ice age (d), pond circularity (e), and center depth together with drainage stage (f). Histograms on the right side show the distribution function of the residuals (b) of all data, (d) separated by the ice age, and (f) separated by drainage stage. Pond circularity is the ratio between the max. PIA disk area and pond area, indicating the shape complexity of the pond.

6.1.2 Accuracy of the in situ data and suggestions for future field campaigns

The accuracy of photogrammetrically derived pond depths was obtained compared to in situ measurements as ground truth. As discussed in König et al. (2020), those measurements also contain substantial sources for errors, namely metersticks that slip into cavities and the complex retracing of in situ measurement points in the airborne data. To compensate for the latter, we averaged above a circular area with a radius of 0.3 m around the in situ points. We strongly recommend a sophisticated system for future campaigns to allow for a direct georeferencing link between in situ and airborne data. A pure GPS-based method proved to be insufficient for constantly drifting ice. Therefore, we recommend a system consisting of GPS base stations with regular freeboard measurements that are recognizable in images and record the geographic position in the Earth-fixed system. Such stations act as optical and geospatial GCPs, also improving the photogrammetric analysis through accurate horizontal and vertical position reference. These reference stations must encompass the measurement area so that in situ measurements within the study area can be efficiently and accurately triangulated between these points. The best choice of tools, e.g., 360° cameras, theodolites, or low-range indoor positioning systems like Bluetooth beacons, needs to be tested. To reduce the error from meterstick point measurements caused by rippled or porous ice, we recommend the usage of a flat plate at the bottom of the depth gauge in future campaigns. This could, however, possibly lead to photogrammetric measurements slightly overestimating the pond depth since the optical transition layer is probably located somewhat deeper than the upper haptic transition between pond water and ice. This requires further testing. It has also been shown on MOSAiC that measurements with echo sounders from remote-controlled boats can be used as a reference. Yet, even acoustic systems, especially in the high-frequency range of 500 kHz as implemented in Böötle on MOSAiC, do not always detect the exact transition from water to ice as shown under laboratory test conditions (Werner, 2023). In general, we assume that in situ reference depths from both manual and acoustic measurements are slightly overestimated by a few centimeters, if at all, and are therefore a suitable means of testing the algorithm for simultaneous flights. The acoustic measurements in this study, however, could only confirm the great depths in the observed ponds and derive trends, but no exact error could be derived due to the time lag of 2 d.

Figure 19 Box-and-whisker plots showing the mean (orange line), lower and upper quartile (box), and the minimum and maximum range (whiskers, 99.3 % of all data points inside this range) of bootstrap samples of the form factor κ that we introduce to upscale single point measurements of pond depth to mean pond depth. Statistics are retrieved from 100 random samples of pond sample size n from the entire set of 1621 ponds used in this study. The shaded area shows the targeted accuracy of the pond form factor as described in the text.

6.1.3 Required flight pattern

Our results suggest that pond depth studies require flight grids with great overlap of adjacent images. The more lateral offset is achieved in addition to forward overlap in the flight direction, the higher the accuracy of the pond depth measurements. However, simultaneously, we showed that small measuring angles, i.e., large aperture angles, should be avoided, as they lead to a greater underestimation of depths due to light refraction. This effect is most likely responsible for the observed discrepancies between the in situ (echo sounder) and aerial pond depths in the deep parts of Mystery Lake on MOSAiC, besides ongoing melting between the observing times. To achieve a high overlap and lateral offset at small measuring angles, aerial images must be captured with a high measurement frequency or low flight speed. To give an example for flight planning, the system we used on MOSAiC is limited to a measurement frequency of 0.25 Hz. To obtain a reasonable ground sampling distance of 5 cm, the camera resolution limits flight altitude to 100 m (300 ft). A circle with the maximum angle of incidence of the measurement of 40° that we defined to prevent errors in the photogrammetric reconstruction from light refraction then has a ground radius of 84 m. It is commonly recommended for photogrammetric measurements to use 80 % forward and 60 % lateral overlap in images for optimal reconstruction. Maximum flight speed would thus be limited to 4.2 m s⁻¹ (8 kts) and the lateral offset of flight lines to 33.6 m. Fortunately, state-of-the-art aerial imaging systems can provide sufficiently higher temporal measurement frequencies to facilitate measurements.

In the photogrammetrically reconstructed topography, the DEM, large-scale gradients, also known as the doming effect, can typically occur. Since the same problem occurred with the PASCAL data, we have introduced tools that include a separate analysis of the ponds so that a derivation of the pond depth is still possible. The effect can be minimized by improving the calibration of the camera model through improved flight patterns, a slightly oblique camera perspective (Wackrow and Chandler, 2008), and, as we found when comparing the PASCAL and MOSAiC data, by not using an additional camera protection window in front of the lens.

6.2 MOSAiC melt pond bathymetry

The second part of the study applies the developed method to data from the MOSAiC campaign. Here, we discuss the major results from that unprecedented insight into temporal and spatial pond evolution and compare it to conventional measurement methods.

6.2.1 Pond bottom and level evolution

The large-scale and high-resolution observation of ponds on the MOSAiC floe points to several previously unconsidered melt season processes. Firstly, large ponds on level ice in the deformed center of the MOSAiC floe were extraordinarily deep (>2 m) before the drainage event. They were thus as deep as the level ice on which they formed was thick (e.g., Itkin et al., 2023; Thielke et al., 2023; von Albedyll et al., 2022), while the pond bottoms themselves still consisted of 1 m and thicker ice, as anecdotal reports from the field showed. After drainage, pond bottoms, even of previously deep ponds, surfaced, although the pond level was close to sea level before. Both these observations are indications for a flexible ice cover (as shown in Fuchs, 2023c) and contradict the traditional assumption of a rigid ice cover on which ponds are only reduced to their lowest points during the drainage event (e.g., Popović et al., 2020).

Secondly, it was mostly assumed or observed so far that pond surfaces are well above sea level before drainage (e.g., Polashenski et al., 2017; Perovich et al., 2021). This was the case only for a few small ponds in the deformed center of the floe. Only their bottom ice seemed impermeable and rigid enough on MOSAiC to resist the hydrostatic head that forms when ponds are above sea level. Even before the drainage event, the largest pond with a water level >0.2 m higher than the sea level had an exposed water surface of 13.4 m in diameter (disk size around the pole of inaccessibility). A total of 46 ponds on the MOSAiC floe were larger, but their pond level was closer to sea level on 30 June. This indicates that large ponds in particular have little possibility of being far above sea level and that blockage processes reducing permeability (Polashenski et al., 2017) are not able to fully sustain the impermeability of the underlying ice. Instead, we speculate that a balance of water inflow and outflow to ponds is established (lateral, vertical, and melting), allowing a pond level only minimally above sea level. Or, alternatively, especially in large ponds, the bottom ice can bend downwards under the load of accumulated meltwater and thus reduce the surface elevation, while the hydrostatic head, causing the bending, is still higher (as discussed in Fuchs, 2023c). However, more detailed investigations of the mechanical properties of the melting ice are required to break down the individual contributions of the processes to the observed results.

The presented observations of the temporal and spatial pond bathymetry evolution indicate that the pond bottoms can bend downwards under the weight of the accumulated meltwater and, after drainage, when the ice is very permeable, bend upwards due to the buoyancy of the ice. This flexibility must be considered when observing hydrostatic heads in the field and parameterizing the pond development in models.

6.2.2 Representativeness of in situ studies

In situ transect studies are among the most common methods to collect representative in situ data during sea ice physical field campaigns. The entire sea ice column, including snow depth, pond depth, and ice thickness, can be probed on walked transect lines, providing comprehensive data on the properties of the atmosphere, snow, ice, and ocean and transition zones between them. However, Perovich et al. (2003) noted that pond coverage along their transect line on the SHEBA campaign (Uttal et al., 2002) with a peak coverage of 40 % was not representative of the entire SHEBA site, for which they retrieved a comparably smaller and more constant pond coverage of less than 24 % from aerial imaging (Perovich et al., 2002). Also, on MOSAiC, differences were found in the derived pond coverage between the walked transect line and high-resolution satellite imagery covering the entire floe (Webster et al., 2022). Webster and colleagues attributed a comparably smaller pond coverage along the transect line to its location close to the floe edge, where lateral runoff of meltwater is known to reduce pond coverage (e.g., Wright et al., 2020). This may also reduce pond coverage on smaller floes as observed by Divine et al. (2015) in the marginal ice zone.

Collocated transect and airborne measurements are rare, as they are only possible with enormous effort on large-scale field campaigns or require high-resolution satellite imagery in clear-sky conditions, which has become available only in recent years. For this reason, combinations of both, such as those listed previously, are rare, and we know of no study in which both methods have been combined with a focus on comparing them instead of broadening the data source. Owing to the data available from Webster et al. (2022), we could therefore, for the first time, include a systematic comparison of derived geometric pond properties between the two measuring methods. The presented results have shown that in situ transect measurements show slightly greater average pond depths than airborne-derived depths and are not entirely representative of the floe before drainage. Combining in situ transect data with high-resolution aerial imaging therefore yields great potential for strengthening studies on the floe scale. The slightly greater observed depths can be expected from manual in situ measurements due to uneven pond bottoms mentioned above. However, we noticed that mean pond depths are apparently commonly derived by averaging depth profiles along single lines through ponds. We want to point out that using a single line to infer the average pond depth is mostly geometrically incorrect due to the mostly round shape of the ponds and possibly overestimates the mean pond depth. Shallower areas at the pond edges are underrepresented in such an extrapolation process.

Because of the possible overestimation of mean pond depth by in situ studies, we investigated the possibility of extrapolating mean pond depths of a large number of ponds from single in situ measurements at their center point. On average, the introduced form factor κ=0.52 indicated a strong correlation between pond depth at the center and mean pond depth. In the center, ponds were about twice as deep as on average, largely independent of other factors such as shape, depth, ice deformation, or stage of evolution. This result, collected from a sample size of 1621 ponds, could help improve mean pond depth estimations in the field where a complete manual assessment of the pond bathymetry is not feasible due to the limited workload, especially for a great number of individual ponds. Measuring only one depth for each of these ponds at the relatively easily findable spot pole of inaccessibility and estimating the mean pond depth using the pond form factor could massively increase the representativeness of the measurement for the mean pond depth.

6.2.3 Large-scale pond observations including satellite data

Despite its proven broad applicability, the presented aerial image method for reconstructing pond geometry remains spatially and temporally limited and requires upscaling by satellite data for pan-Arctic observations. Our surface classifications show no noticeable difference compared to those with another algorithm (OSSP) applied to high-resolution satellite images of the same area (DigitalGlobe WorldView, Webster et al., 2022). Since we also see no further reason for a systematic difference between these high-resolution image input data, we assume that, apart from the different applicability (sky conditions, spatial and temporal coverage), both high-resolution surface class products compare well. Lower resolved optical satellite sensors require spectral unmixing techniques with temporally and spatially varying accuracy, which is discussed and evaluated in detail in other studies (e.g., Niehaus et al., 2023; Lee et al., 2024).

Regarding pond depths, novel approaches are emerging to upscale from a few manual point measurements on the ice to large-scale observations (Farrell et al., 2020; Herzfeld et al., 2023; Buckley et al., 2023; König et al., 2020). The ICESat-2-based observations offer enormous temporal and spatial coverage but are limited to larger ponds. Our high-resolution data can therefore make an essential contribution to evaluating and optimizing them. Similar to Buckley et al. (2023), we found that wide ponds defined here by their diameter at the pole of inaccessibility (which indicates the direction-independent extent of the pond shapes) contributed significantly to the overall melt pond volume on the floe. These ponds were exceptionally deep in our data, which supports the Buckley et al. (2023) statement that ICESat-2 data possibly overestimate the mean pond depth due to the reduced range in the observable melt pond size distribution. With the drainage event, the diameter and, with that, the direction-independent horizontal expansion of open pond water dropped significantly while the pond areas still remained high (caused by braided pond patterns) (Sect. 5.1). This means that many drained ponds observed from the helicopter potentially fell below the minimum size thresholds of the ICESat-2 algorithms, which on the MOSAiC floe could have resulted in an underestimation of the bulk pond volume after drainage by up to 18.7 % to 68.4 % (DDA minimum width) or 92.3 % (UMD-MPA minimum width), assuming that no pond would have been observed from the satellite in a direction longer than that defined by the diameter. While the MOSAiC floe was certainly somewhat unique with its strongly deformed center part, this highlights how much small ponds (of which even the smallest are also cut off due to our 100-pixel threshold) contribute to the overall pond volume, especially after the drainage event when braided pond patterns reduce widespread, lake-like ponds. Incorporating data from the new aerial-based approach into satellite retrievals offers enormous potential for the upscaling of pond properties that should be fully exploited in the future.

7 Conclusions

We proved that pond bathymetry can be accurately derived from a photogrammetrical reconstruction of the ice surface. Aerial images from a monocular airborne camera in motion are sufficient, provided that strong overlap is given between single images. A simple multiplication of the derived water column depth with the refraction index of water (n=1.335) sufficiently corrects measured pond depths for light bending at the water–air interface. This factor naturally varies depending on the angle of incidence of the measurement. Still, here it has been shown that a constant factor is sufficient to be subsequently used in the complex bundle-block adjustment in Agisoft Metashape, which favors nadir measurements. Incident angles at the surface were restricted to smaller than 40° to avoid horizontal alignment errors becoming greater than the ground sampling distance. Deviating from this limitation possibly contributed to a slight underestimation of large depths in the MOSAiC data.

With the newly developed method, we could reconstruct the evolution of pond coverage, depth, and volume on the MOSAiC leg 4 floe. The pond volume of the entire floe was more constant than found by Webster et al. (2022), who extrapolated it from a transect line around the floe. This difference was mainly caused by very deep ponds with large diameters in the floe center that doubled the mean area-specific pond volume of the floe before drainage when coverage was still slightly lower. Harmonization of the pond cover after drainage increased the representativeness of the transect area and redistributed pond volume to ponds with smaller diameters (5 to 20 m), which then contributed significantly (92.5 %) to the overall pond volume. Furthermore, we could detect clear indicators for a flexible ice cover below the ponds and derive a scaling factor for retrieving mean pond depths from single pond depth measurements in the field. We showed that pond depth and volume strongly depend on the pond diameter at the pole of inaccessibility of the pond shape. We therefore recommend this pond diameter as a valuable measure of pond characteristics.

The study showcases unexploited possibilities of aerial imaging of melt ponds. The developed methods and procedures allow pond bathymetry reconstruction solely from one optical camera deployed on a helicopter or airplane. The resulting availability and exact alignment of optical and morphological data provide a unique database that requires only relatively inexpensive instrumentation and evaluation software to be compiled. Although we used the commercial software Agisoft Metashape for the most complex processing steps, we would like to point out that there are open-source projects like MicMac (Rupnik et al., 2017) available, which could further reduce the cost for such systems. Pre-processing and post-processing are completely based on freely available Python and QGIS packages. The only change to previous campaigns, which eventually allows for the pond depth retrieval, is the appropriate design of the flight tracks. This approach is probably directly transferable to uncrewed aerial vehicles (UAVs). Such a transformation from airplane-based observation to UAVs holds great potential for reduced emissions and an economically more friendly collection of sea ice observations. The data acquired during MOSAiC show that the method enables multidimensional studies tackling questions far beyond the current scope.

Appendix A: Calculations

The following geometric approaches were used to derive the horizontal and vertical mismatch caused by refraction in the photogrammetric depth calculations.

A1 Horizontal mismatch

From Fig. 2, one can geometrically derive a set of equations describing the different sections of X.

$$\begin{array}{}\text{(A1)}& {\displaystyle }{X}_{{\mathit{\beta }}_{\mathrm{1}}}& {\displaystyle }=\mathrm{\Delta }X+X_{{\mathit{\alpha }}_{\mathrm{1}}}& {\displaystyle }=Z_{\mathrm{AP}}\cdot \tan {\mathit{\beta }}_{\mathrm{1}}\text{(A2)}& {\displaystyle }{X}_{{\mathit{\beta }}_{\mathrm{2}}}& {\displaystyle }=X_{{\mathit{\alpha }}_{\mathrm{2}}}-\mathrm{\Delta }X& {\displaystyle }=Z_{\mathrm{AP}}\cdot \tan {\mathit{\beta }}_{\mathrm{2}}\text{(A3)}& {\displaystyle }{X}_{{\mathit{\alpha }}_{\mathrm{1}}}& {\displaystyle }=Z_{\mathrm{TP}}\cdot \tan {\mathit{\alpha }}_{\mathrm{1}}\text{(A4)}& {\displaystyle }{X}_{{\mathit{\alpha }}_{\mathrm{2}}}& {\displaystyle }=Z_{\mathrm{TP}}\cdot \tan {\mathit{\alpha }}_{\mathrm{2}}\end{array}$$

Equating these with Z_AP results in

$$\begin{array}{}\text{(A5)}& \mathrm{\Delta }X={\displaystyle \frac{X_{{\mathit{\alpha }}_{\mathrm{2}}}\cdot \frac{\tan {\mathit{\beta }}_{\mathrm{1}}}{\tan {\mathit{\beta }}_{\mathrm{2}}}-X_{{\mathit{\alpha }}_{\mathrm{1}}}}{\mathrm{1}+\frac{\tan {\mathit{\beta }}_{\mathrm{1}}}{\tan {\mathit{\beta }}_{\mathrm{2}}}}}.\end{array}$$

By substituting $X_{{\mathit{\alpha }}_i}$ and applying Snell's law to replace β by the known angle of emergence α, Eq. (A5) becomes

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{\mathrm{\Delta }X}{Z_{\mathrm{TP}}}}={\displaystyle \frac{\tan \left({\mathit{\alpha }}_{\mathrm{2}}\right)\cdot \frac{\tan \left(\arcsin \left(\frac{\sin {\mathit{\alpha }}_{\mathrm{1}}}{n_{\mathrm{water}}}\right)\right)}{\tan \left(\arcsin \left(\frac{\sin {\mathit{\alpha }}_{\mathrm{2}}}{n_{\mathrm{water}}}\right)\right)}-\tan \left({\mathit{\alpha }}_{\mathrm{1}}\right)}{\mathrm{1}+\frac{\tan \left(\arcsin \left(\frac{\sin {\mathit{\alpha }}_{\mathrm{1}}}{n_{\mathrm{water}}}\right)\right)}{\tan \left(\arcsin \left(\frac{\sin {\mathit{\alpha }}_{\mathrm{2}}}{n_{\mathrm{water}}}\right)\right)}}}.\end{array}$$

For this, we use the simplified image of opposite measuring points and two different angles of emergence on the pond surface.

A2 Vertical mismatch

The deviation in depth between measured and true depth is derived from geometric analysis of Fig. 2:

$$\begin{array}{}\text{(A7)}& {\displaystyle }{X}_{\mathrm{2}}-X_{\mathrm{1}}& {\displaystyle }=Z_{\mathrm{AP}}\tan {\mathit{\beta }}_{\mathrm{2}}+Z_{\mathrm{AP}}\tan {\mathit{\beta }}_{\mathrm{1}},\text{(A8)}& {\displaystyle }{X}_{\mathrm{2}}-X_{\mathrm{1}}& {\displaystyle }=Z_{\mathrm{TP}}\tan {\mathit{\alpha }}_{\mathrm{2}}+Z_{\mathrm{TP}}\tan {\mathit{\alpha }}_{\mathrm{1}}.\end{array}$$

Equating both and including Snell's law leads to

$$\begin{array}{}\text{(A9)}& {\displaystyle }\mathit{\gamma }& {\displaystyle }={\displaystyle \frac{Z_{\mathrm{AP}}}{Z_{\mathrm{TP}}}}={\displaystyle \frac{\tan {\mathit{\alpha }}_{\mathrm{1}}+\tan {\mathit{\alpha }}_{\mathrm{2}}}{\tan {\mathit{\beta }}_{\mathrm{1}}+\tan {\mathit{\beta }}_{\mathrm{2}}}},\text{(A10)}& {\displaystyle }\mathit{\gamma }& {\displaystyle }={\displaystyle \frac{\tan {\mathit{\alpha }}_{\mathrm{1}}+\tan {\mathit{\alpha }}_{\mathrm{2}}}{\tan \left(\arcsin \left(\frac{\sin {\mathit{\alpha }}_{\mathrm{1}}}{n_{\mathrm{water}}}\right)\right)+\tan \left(\arcsin \left(\frac{\sin {\mathit{\alpha }}_{\mathrm{2}}}{n_{\mathrm{water}}}\right)\right)}}.\end{array}$$

We apply a limit value analysis for small incident angles. To do so, we first equate both measuring angles α₁=α₂. This results in a simplified form of Eq. (A10):

$$\begin{array}{}\text{(A11)}& \mathit{\gamma }={\displaystyle \frac{\tan \mathit{\alpha }}{\tan \left(\arcsin \left(\frac{\sin \mathit{\alpha }}{n_{\mathrm{water}}}\right)\right)}}.\end{array}$$

In this form, the equation is undefined at zero since tan(0)=0 and arcsin(0)=0. However, with a couple of rearranging tricks, one gets an evaluable formula form. Inserting

$$\begin{array}{}\text{(A12)}& \tan \left(\arcsin \right(x\left)\right)={\displaystyle \frac{x}{\sqrt{\mathrm{1}-x^{\mathrm{2}}}}}\end{array}$$

into Eq. (A11) leads to

$$\begin{array}{}\text{(A13)}& \mathit{\gamma }={\displaystyle \frac{\tan \left(x\right)}{\left(\frac{\frac{\sin \left(\mathit{\alpha }\right)}{n_{\mathrm{water}}}}{\sqrt{\mathrm{1}-{\left(\frac{\sin \left(\mathit{\alpha }\right)}{n}\right)}^{\mathrm{2}}}}\right)}}.\end{array}$$

By rearranging and the definition of the tangent we obtain

$$\begin{array}{}\text{(A14)}& \mathit{\gamma }\left(\mathit{\alpha }\right)=\cos \left(\mathit{\alpha }\right)\cdot \sqrt{n_{\mathrm{water}}^{\mathrm{2}}-{\sin }^{\mathrm{2}}\left(\mathit{\alpha }\right)}.\end{array}$$

Appendix B: MOSAiC pond properties resolved by pond area

Figure B1 Distribution functions of pond volume (a), pond areal fraction (b), mean pond depth (c), and number of ponds (d) in 21 pond area bins from 0 to 10 000 m² and larger than 10 000 m² on the MOSAiC floe. Panel (e) shows the number of ponds resolved by pond diameter as used in Fig. 15. Shown are data retrieved from photogrammetric reconstructions from aerial images collected on survey flights on 30 June 2020 (pre-drainage), 17 July, and 22 July 2020 (post-drainage).