2.1. Snow layer tracking through deep learning

Feature extraction capabilities and good generalizability of deep learning networks have been used on Snow Radar images to track snow layers in Yari et al. (Reference Yari, Rahnemoonfar, Paden, Oluwanisola, Koenig and Montgomery2019), Yari et al. (Reference Yari, Rahnemoonfar and Paden2020), Rahnemoonfar et al. (Reference Rahnemoonfar, Yari, Paden, Koenig and Ibikunle2021), Varshney, Rahnemoonfar, Yari, Paden, et al. (Reference Varshney, Rahnemoonfar, Yari, Paden, Ibikunle and Li2021), Varshney, Yari, et al. (Reference Varshney, Yari, Chowdhury and Rahnemoonfar2021), and Wang et al. (Reference Wang, Xu, Paden, Koenig, Fox and Crandall2021). Deep learning has also been used to segment the radar depth sounder (RDS) images in Donini et al. (Reference Donini, Bovolo and Bruzzone2022), Ghosh and Bovolo (Reference Ghosh and Bovolo2022a), and Ghosh and Bovolo (Reference Ghosh and Bovolo2022b) into ice and bedrock through various modifications of UNet (Ronneberger et al., Reference Ronneberger, Fischer and Brox2015), with or without transformers (Vaswani et al., Reference Vaswani, Shazeer, Parmar, Uszkoreit, Jones, Gomez, Kaiser, Polosukhin, Guyon, Luxburg, Bengio, Wallach, Fergus, Vishwanathan and Garnett2017). In this paper, we focus on the shallow sensor, the Snow Radar (Gogineni et al., Reference Gogineni, Yan, Gomez, Rodriguez-Morales, Paden and Leuschen2013), which tracks the snow layers.

Wang et al. (Reference Wang, Xu, Paden, Koenig, Fox and Crandall2021) used a CNN-RNN (recurrent neural network) architecture for tiered segmentation of radar images to identify the top 4–6 snow layers. Similarly, Ibikunle et al. (Reference Ibikunle, Paden, Rahnemoonfar, Crandall and Yari2020) built an iterative neural network architecture to track snow layers, but on simulated radar images, as real radar images were too complex for feature extraction. Varshney, Rahnemoonfar, Yari, Paden, et al. (Reference Varshney, Rahnemoonfar, Yari, Paden, Ibikunle and Li2021) used pyramid pooling modules to learn local-to-global spatio-contextual information of snow layer pixels and perform semantic segmentation. In the line of multi-scale architectures, Yari et al. (Reference Yari, Rahnemoonfar, Paden, Oluwanisola, Koenig and Montgomery2019) first explored multi-scale contour detection CNNs for the purpose of snow layer extraction. The authors noted that using popular pre-training strategies would not work on Snow Radar images due to the inherent noise present in these images. Further extensions of this work by Rahnemoonfar et al. (Reference Rahnemoonfar, Yari, Paden, Koenig and Ibikunle2021) and Yari et al. (Reference Yari, Rahnemoonfar and Paden2020) showed the usefulness of multi-scale architectures on synthetic radar images and temporal transfer learning, respectively. Subsequently, wavelet-based multi-scale architectures for snow layer tracking were first developed by Varshney, Yari, et al. (Reference Varshney, Yari, Chowdhury and Rahnemoonfar2021) which showed preliminary results in this domain. This prior work used a VGG-13 architecture (Simonyan & Zisserman, Reference Simonyan and Zisserman2015) on a small Snow Radar dataset, which achieved an F-score of a little over 0.7. We expand upon this work by training a VGG-16 architecture, which is known to give robust representations (Rahnemoonfar et al., Reference Rahnemoonfar, Yari, Paden, Koenig and Ibikunle2021), on a new dataset specifically catered towards dry zone areas, where the snow layers can be tracked. Furthermore, and most importantly, we also calculate the individual layer depths in terms of the number of pixels. We also show that using wavelet transforms of each scale of a multi-scale network helps in learning inherent image features, as compared to using multi-level wavelet transforms of the input image.

2.2. CNNs with wavelet transformations

Wavelet transforms have had immense applications in image denoising (Bnou et al., Reference Bnou, Raghay and Hakim2020; P. Liu et al., Reference Liu, Zhang, Lian and Zuo2019; Mallat, Reference Mallat1989), image compression (Naveen Kumar et al., Reference Naveen Kumar, Jagadale and Bhat2019; Xiong & Ramchandran, Reference Xiong, Ramchandran and Bovik2009), and image restoration (Bae et al., Reference Bae, Yoo and Ye2017; Huang et al., Reference Huang, He, Sun and Tan2017; P. Liu et al., Reference Liu, Cheng, Hu, Bian, Zhang, Bai and Tang2019). Many recent works such as P. Liu et al. (Reference Liu, Zhang, Lian and Zuo2019), Williams and Li (Reference Williams and Li2018), Huang et al. (Reference Huang, He, Sun and Tan2017), Han and Ye (Reference Han and Ye2018), and Bae et al. (Reference Bae, Yoo and Ye2017) have used the downsampling properties of wavelets and replaced the convolutional or pooling layers in a CNN, since incorporating the wavelet coefficients typically helps reduce information loss (P. Liu et al., Reference Liu, Cheng, Hu, Bian, Zhang, Bai and Tang2019). Williams and Li (Reference Williams and Li2018) showed that using wavelet transform for downsampling would create much cleaner and sharper images, as compared to using pooling layers, and improve generalizability. Similarly, P. Liu et al. (Reference Liu, Zhang, Lian and Zuo2019) found that wavelet transforms not only enlarge a kernel’s receptive field but also prevent information leakage which generally takes place during a pooling operation. By using multi-level wavelet transforms, Han and Ye (Reference Han and Ye2018) developed deep convolutional framelets to reconstruct sparse-view computed tomography images and Huang et al. (Reference Huang, He, Sun and Tan2017) developed a wavelet-based loss function to super-resolve facial images. Furthermore, Bae et al. (Reference Bae, Yoo and Ye2017) showed that residual learning improves by training on wavelet subbands. The variety of these works shows that by incorporating wavelet transforms in a CNN, the feature extraction capability of the latter improves due to sharper feature maps, increased receptive field, and enhanced residual learning.

3.1. Wavelet transform

We use the discrete wavelet transform (DWT) for our computations. For an image $ f\left(x,y\right) $ having dimensions $ X\times Y $ , DWT is given as Equations 3.1 and 3.2 for $ j\ge {j}_0 $ :

(3.1) $$ {W}_{\phi}\left({j}_0,m,n\right)=\frac{1}{\sqrt{X.Y}}\sum \limits_{x=0}^{X-1}\sum \limits_{y=0}^{Y-1}f\left(x,y\right){\phi}_{j_0,m,n}\left(x,y\right) $$

(3.2) $$ {W}_{\psi}^i\left(j,m,n\right)=\frac{1}{\sqrt{X.Y}}\sum \limits_{x=0}^{X-1}\sum \limits_{y=0}^{Y-1}f\left(x,y\right){\psi}_{j,m,n}^i\left(x,y\right) $$

In these equations, $ {W}_{\phi } $ is the approximation coefficient (A), and $ {W}_{\psi}^i $ are the detail coefficients for each level of wavelet transform $ j $ , where $ i\in \left\{H,V,D\right\} $ . H, V, and D are the horizontal, vertical, and diagonal details, respectively, $ m,n $ are the subband dimensions (Williams, Li, et al., Reference Williams and Li2018), and $ {j}_0 $ is an arbitrary starting level. Each subsequent level of wavelet transform is computed on the previous level’s approximation coefficient. In the beginning, the input raw image is treated as an approximation coefficient. Furthermore, $ \phi $ is a scaling function, and $ \psi $ is the wavelet function, both of which downsample an input image by a factor of 2 in both X and Y dimensions. Readers are encouraged to go through Daubechies (Reference Daubechies1992) to get a detailed understanding of $ \phi $ , $ \psi $ , and the wavelet basis functions. We use the downsampling property of wavelets and combine their transforms with the side outputs of multi-scale networks. Figure 4 shows a level 2 (i.e. j = 2) transform of an input image. In our work, we experiment with three popular discrete wavelets, i.e. Haar, Daubechies, and “Discrete” Meyer. We will be abbreviating them as “haar,” “db,” and “dmey” in subsequent sections for ease of reading. The wavelet functions ( $ \psi $ ) for the three wavelet types have been plotted in Figure 5.

Figure 4. A level 2 wavelet transform of a given input image. The subscript denotes the level number.

Figure 5. The three types of wavelets that we use for our experiments, and their shorthand representation in parenthesis.

3.2. Multi-scale CNN (MS-CNN)

In this work, we use a multi-scale network built for contour detection (Xie & Tu, Reference Xie and Tu2015), which was proven to be useful for snow layer tracking (Rahnemoonfar et al., Reference Rahnemoonfar, Yari, Paden, Koenig and Ibikunle2021; Yari et al., Reference Yari, Rahnemoonfar, Paden, Oluwanisola, Koenig and Montgomery2019). This is a VGG-16 (Simonyan & Zisserman, Reference Simonyan and Zisserman2015) network with the terminal fully connected layers and the last pooling layer removed. The final convolutional layers of every stage, right before the max pooling layer of every stage, are convolved with a 1 × 1 filter to generate what we call “side output(s)”. Each side output is at a downsampled resolution of the previous stage’s side output. All side outputs are finally upsampled through transposed convolutions, followed by cropping, to match the resolution of the input image. The five side outputs, from the five stages, are concatenated together to form a “fuse” layer. All five side outputs along with the final “fuse” layer are trained together in a deeply supervised manner (Lee et al., Reference Lee, Xie, Gallagher, Zhang, Tu, Lebanon and Vishwanathan2015) with a cumulative loss function, Equation 3.5, as explained in the following paragraph.

To train the network, we calculate the binary cross entropy loss ( $ l $ ) for every pixel $ i $ as follows:

(3.3) $$ l\left({x}_i;W\right)=\left(\begin{array}{ll}\alpha \cdot \log \left(1-{x}_i\right)& if\;{y}_i=0\\ {}\beta \cdot \log \left({x}_i\right)& if\;{y}_i=1\end{array}\right. $$

In Equation 3.3, $ {x}_i $ is the sigmoid activation map obtained from a network with weights $ W $ , $ {y}_i $ is the ground truth label of the corresponding pixel in the input image $ I $ having a total of $ \mid I\mid $ pixels, and $ \alpha $ , $ \beta $ are defined as Equation 3.4.

(3.4) $$ {\displaystyle \begin{array}{c}\alpha =\lambda \cdot \frac{\mid {Y}^{+}\mid }{\mid {Y}^{+}\mid +\mid {Y}^{-}\mid}\\ {}\beta =\frac{\mid {Y}^{-}\mid }{\mid {Y}^{+}\mid +\mid {Y}^{-}\mid}\end{array}} $$

In Equation 3.4, $ \mid {Y}^{+}\mid $ denotes the count of all positive labels, i.e. those pixels representing the top of a layer ( $ {y}_i=1 $ ) and $ \mid {Y}^{-}\mid $ denotes the count of all negative labels i.e. all other pixels which are background ( $ {y}_i=0 $ ). $ \lambda $ is a hyperparameter used to balance these positive and negative labels.

The total loss is computed as Equation 5 where $ k $ depicts each of the side outputs or stages, i.e. $ K=5 $ . This network, called Multi-Scale CNN, is the current state-of-the-art model for snow layer tracking (Rahnemoonfar et al., Reference Rahnemoonfar, Yari, Paden, Koenig and Ibikunle2021) and forms our baseline architecture to compare with our wavelet-based architectures. Multi-Scale CNN is shown in Figure 6 and most of the time abbreviated as MS-CNN for the rest of the paper.

(3.5) $$ L(W)=\sum \limits_{i=1}^{\mid I\mid}\left(\sum \limits_{k=1}^Kl\left({x}_i^k;W\right)+l\Big({x}_i^{fuse};W\Big)\right) $$

Figure 6. MS-CNN—The multi-scale architecture of Rahnemoonfar et al. (Reference Rahnemoonfar, Yari, Paden, Koenig and Ibikunle2021) which forms our base model.

3.3. Wavelet combined neural networks

We set up two wavelet-ased architectures, and exploit the multi-level nature of wavelet transforms to embed them into the base multi-scale architecture, MS-CNN. The two wavelet-based network modifications are termed as “WaveNet” and “Skip-WaveNet,” and are described below:

3.3.1. WaveNet

In Figure 3, from stages 2 to 5 of the base architecture, the feature maps get downsampled to a scale of 2× to 16× (with respect to the input image) to form the side outputs. Multi-level wavelet transforms of the input image also give us downsampled features at similar scales, which are sharp, denoised, and contain the horizontal, vertical, and diagonal detail coefficients (see subsection 3.1). By fusing wavelet information of the image to every scale, local as well as global, of the multi-scale network, the side outputs can be enriched for feature extraction. Hence, in this architecture, we take a level 4 wavelet transform (since there are four stages which have downsampled resolution as compared to the input image—stages 2 to 5) of the input image and fuse the detail coefficients in the following manner: we concatenate all three detail coefficients i.e. $ H,V,D $ of a wavelet transform of level $ l $ of the input radar echogram to the side output $ l+1 $ of the base architecture. This means that the detail coefficients from level 1 of the wavelet transform fuse with side output 2 of the architecture, detail coefficients from level 2 of the wavelet transform fuse with side output 3 of the architecture, and so on. The WaveNet architecture is shown in Figure 7.

Figure 7. WaveNet—A wavelet-based architecture. Here, the input image goes through a multi-level wavelet transform, where each level is shown with a different color, ranging from light pink to yellow.

3.3.2. Skip-WaveNet

The side outputs of MS-CNN are feature maps containing image patterns extracted at different scales, from local to global. These feature maps contain an improved representation of the features, i.e. are “more learned,” after every iteration of training and also contain noise in the form of unuseful network weights. Hence, they have the potential to be denoised further which can add value to overall network training. So for the Skip-WaveNet architecture, we take a wavelet transform of each side output and add it as an extra layer to the successive side output to supplement network learning. This will help form skip connections, support residual learning (Bae et al., Reference Bae, Yoo and Ye2017; He et al., Reference He, Zhang, Ren and Sun2016) and also propagate information between scales. Since the wavelet transform is at a downsampled resolution to its input signal, the dimensions of a transform of side output $ s $ will be the same as the dimensions of side output $ s+1 $ . Furthermore, by fusing the wavelet transform of a scale to its successive scale, we are also propagating information between scales which can support network learning. With this intuition, we form the Skip-WaveNet architecture as follows: we take a level 1 wavelet transform of a side output $ s $ ( $ s\in \left[1,4\right] $ ) and fuse its detail coefficients to side output $ s+1 $ . In this case, for a given input image, the wavelet transforms will be renewed after each epoch of training. This can be compared against the WaveNet architecture where the transforms generated for an input image will be the same during all epochs. Hence, such multi-scale “dynamic” wavelets should improve the denoising and feature extraction capabilities of the network. The architecture of Skip-WaveNet is shown in Figure 8.

Figure 8. Skip-WaveNet—A wavelet-based architecture with skip connections. Here, we take a level 1 wavelet transform of each side output. This is in contrast to WaveNet where the wavelet transform was of the input image.

Apart from the above-defined modifications, WaveNet and Skip-WaveNet are trained in the same way as the base architecture MS-CNN, i.e. they are trained with deep supervision through binary cross-entropy loss, and have the final fusion of five side outputs, etc.

3.4. Experimental setup

All the architectures, including the base architecture, and wavelet combined neural networks share the same hyperparameters as those used in Xie and Tu (Reference Xie and Tu2015) and Rahnemoonfar et al. (Reference Rahnemoonfar, Yari, Paden, Koenig and Ibikunle2021). For our experiments, we augment the training dataset with scaling factors $ \in \left[\mathrm{0.25,0.5,0.75}\right] $ and left-to-right flipping. The augmentation helps in expanding the training dataset to 6,430 images, which we train for 15 epochs. All networks took approximately 10 hours, on average, when trained on NVIDIA GeForce RTX 2080 Ti GPU with an Intel Core i9 processor.

3.5. Post-processing and evaluation

We perform non-maximum suppression (NMS) on the network outputs to make finer predictions of the layers. NMS gives grayscale predictions, which we compare against our ground truth by evaluating the ODS (optimal dataset scale) and OIS (optimal image scale) F-scores (Y. Liu et al., Reference Liu, Zhang, Lian and Zuo2019). The ODS F-score looks for a single threshold across the entire dataset which can binarize the predicted image, and give the most optimum F-score with respect to the ground truth. On the other hand, the OIS F-score searches for a similar optimum threshold for each image. The optimum F-scores thus found from each image are then averaged over the entire dataset to obtain the OIS F-score. Here, the F-score refers to the standard F1-score used in computer vision, defined as Equation 3.6. In this equation, TP, FP, and FN are true positives, false positives, and false negatives, respectively. We also use average precision (AP), which is the area under a precision-recall curve (Su et al., Reference Su, Yuan and Zhu2015), as an evaluation metric of network performance.

(3.6) $$ F=\frac{TP}{TP+\frac{1}{2}\left( FP+ FN\right)} $$

Based on the optimum threshold we get, we binarize the network outputs to get predictions similar to the ground truth. We then calculate the row indices (or the “range bin” indices of an echogram) of each layer in the predicted image matrix and compare them with corresponding row indices of the manually labeled ground truth to evaluate the mean absolute error (MAE) of layer depth using Equation 3.7. The overall MAE was then obtained by estimating the MAE for each layer in every echogram for the entire test set using Equation 3.8.

(3.7) $$ {MAE}_j^k=\frac{\sum_{i=1}^{i={W}_k}\mid {GTD}_i^{k,j}-{PD}_i^{k,j}\mid }{W_k} $$

(3.8) $$ MAE=\frac{1}{T}\sum \limits_{k=1}^{k=T}\frac{\sum_{j=1}^{j={N}_k}{MAE}_j^k}{N_k} $$

In Equations 3.7 and 3.8, $ GTD $ , and $ PD $ are the ground truth and predicted layer depths, respectively, in terms of the number of pixels. For example, $ {PD}_i^{k,j} $ represents the the row index of the $ {j}^{th} $ predicted layer in the $ {i}^{th} $ column of the echogram $ k $ . $ {W}_k $ and $ {N}_k $ _, respectively, represent the total number of columns (width), and the total number of traced layers, in an echogram $ k $ . $ T $ represents the total number of echograms in the test set.

3.6. Comparison with in-situ stake measurements

The snow thickness and annual accumulation rates can be computed from traced annual layers. To assess the effects of tracing imperfections on the annual accumulation rate estimates, we compare the radar-derived accumulation rates with the ground measurements. Of the available in-situ stake measurements compiled by the SUMup group (see Section 4), the closest measurement to the flight line of the radar data set used in this paper is chosen as the reference for comparison. The stake location is only $ \sim $ 16 km away and is thus suitable for reliable comparison. Given the scale of the ice-sheet, we do not expect the accumulation to vary at such a short distance. We estimate the radar-derived snow accumulation rates similar to J. Li et al. (Reference Li, Rodriguez-Morales, Fettweis, Ibikunle, Leuschen, Paden, Gomez-Garcia and Arnold2023) which uses Benson’s density-depth model (Benson, Reference Benson1960) and an interpretation model from Clarke et al. (Reference Clarke, Cross and Benson1989). The in-situ stake measurements provide the required initial values of snow density and accumulation at the surface for the interpretation model (Clarke et al., Reference Clarke, Cross and Benson1989). The monthly stake-measured snow accumulations between September and August in two consecutive years are added up as the annual accumulation rate. The MAE effects are assessed by comparing the radar-derived accumulation rates with and without the MAE of each layer, which is calculated by Equation 3.8.

5.1. Wavelet combined neural networks vs MS-CNN

We tabulate the ODS and OIS F-scores, the average precision (AP), as well as the mean absolute error (MAE) of depth estimates in Table 2. These estimates are calculated on the final fuse layer of our experiments. From Table 2, we see that all wavelet combined neural networks, except for WaveNet with a dmey wavelet, perform better than the baseline MS-CNN model. The wavelet-based networks also give higher AP and MAE scores (3.31 pixels) as compared to CNN3B+RNN (Wang et al., Reference Wang, Xu, Paden, Koenig, Fox and Crandall2021). Wang et al. (Reference Wang, Xu, Paden, Koenig, Fox and Crandall2021) were able to predict the top four to six layers from an echogram, while the multi-scale architecture, whether MSCNN or wavelet-based, can trace out the deeper, twenty to thirty, layers from the echogram (Figures 1-3).

Table 2. Evaluation metrics were obtained across different model architectures on the test set. Highest scores are highlighted in bold

In the wavelet-based networks, Skip-WaveNet with a dmey wavelet achieves the highest F-scores, with all Skip-WaveNet models performing better than WaveNet. In the rest of the manuscript, we use the words Skip-WaveNet and Skip-WaveNet-dmey interchangeably. The wavelet transforms help in improving the ODS F-score by 2.11%–3.99% and the OIS F-score by 1.96%–3.7%, over the baseline architecture, for different experiments. Thus, the wavelet transforms not only help in denoising feature maps but also extract features at the global scale (fifth stage), which the comparative MS-CNN does not. Since images are extremely noisy with different images having varying levels of intensity and contrast (see the different sample radar echograms in Figures 1-3), the layers can be almost impossible to see by human eyes. Hence, extracting features at a global fifth stage is a useful property that the wavelet architectures have achieved. In Figures 1-3, it can be seen that the layers detected by the wavelet network give a consistent and useful tracking result than the baseline MS-CNN model. The baseline model is misguided by a lot of degraded features in the input radar echograms, and falsely classifies redundant backscatterings as a layer, as can be seen by its intermediary predictions between two layers, highlighted in red in the figures. The layers predicted by MS-CNN mostly have missing values and are discontinuous in nature (i.e. have gaps in between) giving a reduced tracking accuracy (Table 2). MSCNN’s MAE is worse by $ \sim $ 6 pixels which will give inaccurate thickness estimates eventually adding to uncertainties in accumulation rates (Kahle et al., Reference Kahle, Steig, Jones, Fudge, Koutnik, Morris, Vaughn, Schauer, Stevens and Conway2021). Further, it should be noted that all multiscale architectures (MS-CNN, WaveNet, and Skip-WaveNet) have the same baseline CNN architecture, and hence the same number of trainable parameters. In WaveNet and Skip-WaveNet, the wavelet coefficients are additional filters computed from either the convolutional layers (side outputs) or the input images, and hence neither act as parameters which require training nor are updated by backpropagation. In this regard, both Skip-WaveNet and WaveNet are able to achieve a higher tracking accuracy even though their complexity, in terms of number of trainable parameters, is the same as that of MS-CNN.

5.2. Skip-WaveNet vs WaveNet

Our experiments confirm that wavelet-based architectures, by supplying additional information from the wavelet coefficients, enhance feature extraction. From Table 2 we see that Skip-WaveNet performs better than others, where the wavelet transform obtained in each epoch is from a learned and improved side output (so the wavelets are “dynamic” and change after each epoch), whereas in WaveNet, the multi-level wavelet transforms are always fixed for a given input image for every epoch (the wavelets are in a way “static” across the entire training period). We further see that the performance of the WaveNet experiment deteriorates with the wavelet complexity, whereas the performance of Skip-WaveNet improves with wavelet complexity (ODS F-score of dmey $ > $ db $ > $ haar). What would be interesting for future studies is to investigate how the different types of wavelets are affecting network learning in different ways.

5.3. Validation with snow stake measurements

In addition to the geographical proximity, the measurements at the selected snow stake site (Dibb & Fahnestock, Reference Dibb and Fahnestock2004; the cyan dot in the left panel of Figure 9) span the temporal range of 2003 to 2016, thereby exhibiting a temporal overlap with the radar data collection period. This temporal alignment renders it particularly well-suited for conducting comparative analyses with our model predictions applied to the 2012 IceBridge dataset. We use this snow stake site and find the test echogram closest to it, to compare Skip-WaveNet’s predictions over this echogram and the measurements at the selected stake site. Figure 11 shows the annual accumulation rate for different years (2004–2012), as well as the associated error in these measurements at the stake measurement site. Similarly, the plot also shows the radar-derived accumulation rate and its error due to the mean absolute error achieved by SkipWaveNet-dmey for the first eight annual snow layers between 2004 and 2011. The depth of the 2004 layer is ~4.5 m below the surface. As can be seen from the plot, the estimates from SkipWaveNet are within the error bounds of the stake accumulation measurements. Moreover, the radar-derived accumulation rates from training labels (i.e. without MAE, and marked by *) agree well with the accumulation rate measured at the stake site. The MAE of the first eight layers in SkipWaveNet’s predictions is 2.2 pixels corresponding to an error in accumulation rate estimation of $ \sim $ 0.011 m w.e. a-1. Therefore the tracing improvement of 6.2 pixels by Skip-WaveNet dmey over MS-CNN (see Table 2) results in reducing the uncertainty by 0.031 m w.e. a-1 in accumulation rate estimation. We also compare the radar-derived accumulation rates from the mean absolute error of SkipWaveNet and MSCNN for the same echogram in Figure 12 which shows that Skip-WaveNet was able to give a smaller error in accumulation rates estimates as compared to MS-CNN.

Figure 11. Comparison of the accumulation rate and its error at the site of stake measurements (blue) with the radar-derived accumulation rate from the traced annual layers marked by the training labels (*), and with the mean absolute error achieved by SkipWaveNet (red).

Figure 12. Comparison of the radar-derived accumulation rate from the mean absolute error of SkipWaveNet (red) and MSCNN (black), both centered at the radar-derived accumulation rate from the traced annual layers (*). Clearly, SkipWaveNet gives a lower error.

At this stake site, however, the available density measurements were made between 1997 and 2002 in different months. In Figure 13 we plot these density measurements, made from the surface to a depth of 3 m, in the form of blue circles. The red line shows the density-depth profile obtained from linear fitting with the least mean square error. The orange line shows the density-depth profile used by the interpretation model (Clarke et al., Reference Clarke, Cross and Benson1989) for the accumulation rate estimation from tracked snow layers in a radar echogram. Figure 14 presents the corresponding profiles of snow deposition age and dielectric constant versus depth, respectively.

Figure 13. Snow density measurements at the stake site (Dibb & Fahnestock, Reference Dibb and Fahnestock2004) in blue, their linear fit in red, and the density-depth profile (through the interpretation model of Clarke et al. (Reference Clarke, Cross and Benson1989)) of snow layers traced from the closest radar echogram in orange.

Figure 14. Profiles of snow layer deposition age (in blue) and dielectric constant (in red) versus depth estimated through the interpretation model of Clark e et al. (1989) at the location of the radar echogram closest to the stake site (Dibb & Fahnestock, Reference Dibb and Fahnestock2004). Layers traced by Skip-WaveNet on this radar echogram are marked by circles.