Rimita Lahiri1, Tiantian Feng1, Rajat Hebbar1, Catherine Lord2, So Hyun Kim3,
Shrikanth Narayanan1

# Robust Self Supervised Speech Embeddings for Child-Adult Classification in
Interactions involving Children with Autism

###### Abstract

We address the problem of detecting who spoke when in child-inclusive spoken
interactions i.e., automatic child-adult speaker classification. Interactions
involving children are richly heterogeneous due to developmental differences.
The presence of neurodiversity e.g., due to Autism, contributes additional
variability. We investigate the impact of additional pre-training with more
unlabelled child speech on the child-adult classification performance. We pre-
train our model with child-inclusive interactions, following two recent self-
supervision algorithms, Wav2vec 2.0 and WavLM, with a contrastive loss
objective. We report $9-13\%$ relative improvement over the state-of-the-art
baseline with regards to classification F1 scores on two clinical interaction
datasets involving children with Autism. We also analyze the impact of pre-
training under different conditions by evaluating our model on interactions
involving different subgroups of children based on various demographic
factors.

Index Terms: speech, child-adult classification, self-supervision, autism

##  1 Introduction

Autism Spectrum Disorder (ASD) is a neuro-developmental disorder,
characterized by deficits in social and communicative abilities along with
restrictive repetitive behavior [1, 2]. Individuals with ASD tend to show
symptoms of anomalies in language, non-verbal comprehension, expressions and
vocal prosody patterns [3, 4]. In the United States, the prevalence of ASD in
children has steadily increased from 1 in 150 [5] in 2002 to 1 in 44 in 2022.
It is critical to develop early ASD diagnosis to create timely interventions.
One of the most common observation tools supporting ASD diagnostic and
intervention efforts includes clinically-administered semi-structured dyadic
interactions between the child and a trained clinician [6, 7]. Computational
analysis of such interactions provides evidence-driven opportunities for the
support of behavioral stratification as well as diagnosis and personalized
treatment.

Figure 1: Schematic overview of the proposed two-step recipe for child-adult
classification

However, with regards to behavioral feature extraction and analysis for these
dyadic interactions, prior works have primarily relied on human-annotated data
segmentation by speaker labels, which is expensive and time-consuming to
obtain, especially for large corpora. Computational modeling of naturalistic
conversations has gained a lot of attention in the past few decades because of
its potential in rich human behavioral phenotyping. Hence, it is desirable to
conduct automatic analysis of these interactions using signal processing and
machine learning. Specifically, one fundamental module for supporting
automated processing of child-adult interactions is the task of child-adult
speech classification i.e., distinguishing the speech regions of the child
from those of an interacting adult. Analysis of child speech is more
challenging than adult speech because of the wide variability and
idiosyncrasies associated with child [8, 9, 10]. An additional layer of
complexity arises while analyzing speech for the clinical domain, as different
clinical conditions may lead to unique patterns in language and speech, making
it challenging for current computational approaches to capture.

Training a robust child-adult classifier is challenging for two main factors:
the scarcity of reliably labeled datasets containing child speech and the
larger within-class variability due to the changes in child speech based on
demographic factors like age, gender, and developmental status including any
clinical symptom severity [11]. Most recent works addressing the problem of
speaker diarization have primarily targeted fine-tuning the pre-trained models
by optimizing a supervised objective. So far, Self-Supervised Learning (SSL)
algorithms are largely under-explored for leveraging unlabelled child speech
for developing speaker discriminative embeddings, especially in real-world
settings such as clinical diagnostic and monitoring sessions. Specifically,
there is a limited understanding of how the performance of these models varies
across children with different demographics, including age and gender.

Contributions of this paper: We address the above questions by evaluating the
impact of including more child speech, during pre-training on the child-adult
speaker classification. We choose Wav2vec 2.0 (W2V2)-base and WavLM-base+ as
the backbone models. The detailed contributions of this work are summarized
as:

  * •

Our work represents one of the first attempts to leverage unlabelled child
speech in pre-training for developing speaker discriminative embeddings,
especially due to the vast inherent heterogeneity in the data arising from
developmental differences.

  * •

We experimentally substantiate the effectiveness of our method for downstream
child-adult speaker classification tasks using W2V2 and WavLM and report over
13% and 9% relative improvement over the base models in terms of F1 scores in
two datasets, respectively.

  * •

We also illustrate and analyze the performance of the proposed method among
different subgroups of children based on demographic factors.

##  2 Background

###  2.1 Self-supervision in speech

The need for building speech processing frameworks in low/limited resource
scenarios has spurred significant efforts on unsupervised, semi-supervised and
weakly supervised learning strategies to reduce reliance on labeled datasets.
The success of SSL [12] in natural language processing, notably due to its
generalizability and transferability, has also inspired its adoption within
the speech domain. Early studies explored SSL in speech with generative loss
[13, 14], while more recent ones have focused on discriminative loss [15, 16]
and multi-task learning objectives [17, 18]. The current approach in this
realm follows a two-step process: first pre-train a model in a self-supervised
manner on large amounts of unlabeled data to encode general-purpose knowledge,
and next specialize the model on various downstream tasks through fine-tuning.
Past studies have reported the efficacy of SSL algorithms by leveraging the
pre-trained embeddings on downstream tasks including ASR [15], speaker
verification [19], speaker identification [20], phoneme classification [21],
emotion recognition [22], spoken language understanding [22], and TTS [23].

###  2.2 Child-adult classification in the ASD domain

Child-adult classification is among the more difficult tasks within speaker
diarization, due to the challenges related to "in the wild" child speech in
naturalistic conversational settings including short speaker turns, varied
noise sources and a larger fraction of overlapping speech. Early diarization
solutions involving child speech used traditional feature representations
(MFCCs, PLPs) [24, 25]. In [26], the authors introduced several methods for
processing audio collected from children with autism using a wearable device.
Later, deep speech representations, i-vectors [26] and x-vectors [27] were
studied for this task. A variety of challenges, both from signal processing
and limited data availability, have been identified and addressed. In [11],
the authors have proposed an adversarial training strategy to address the
large within- and across-age and gender variability due to developmental
changes in children. Alternatively, in [28], pre-trained x-vectors were fine-
tuned for child/adult speaker diarization using a meta-learning paradigm,
namely prototypical networks. Moreover, the role of the amount of child speech
in building deep neural speaker representations was studied in [29] and their
experimental results confirm that including more child data indeed enhances
the task performance in a supervised setup.

##  3 Datasets

Our child-inclusive data come from interactions in a clinical setting,
specifically obtained during the administration of two clinical protocols
related to developmental disorders. The first protocol is the gold standard
Autism Diagnostic Observation Schedule (ADOS) [6], used for diagnostic
purposes. The second protocol is a recently proposed outcome-measure focused
instrument Brief Observation of Social Communication Change (BOSCC) [7] for
tracking changes in social and communicative skills during the course of
treatment. A typical ADOS session lasts $40-60$ minutes and contains multiple
(usually $10-15$) semi-structured activities for addressing specific symptoms
related to ASD. Usually these interactions aim to elicit spontaneous responses
from children under different circumstances to obtain a diagnostic score for
classifying children with and without ASD. A BOSCC session is usually 12
minutes long, consisting of two $2min$ conversational talk sessions and two
$4min$ play sessions where the child plays with a toy.

In our pre-training experiments, we use a dataset consisting of 369 recordings
of unlabelled BOSCC sessions comprising approximately $100K$ utterances. For
the fine-tuning experiments, we use two different corpora, ADOSMod3 and
Simons. The ADOSMod3 corpus was collected across 2 clinical sites. These data
are from administrations of the ADOS Module-3 designed for verbally fluent
children, with a focus on the Social Difficulties and Annoyance and Emotional
sub-tasks for this work. The data consist of total 346 sessions collected from
$165$ children ($86$ ASD, $79$ Non-ASD). The Simons corpus used in our study
consists of a combination of clinically administered ADOS ($n=6$) and BOSCC
($n=33$) sessions collected across 4 sites and these sessions were labeled by
trained annotators to extract speaker timestamps. The details of datasets are
reported in Table 1.

Table 1: Session-level statistics of child-adult corpora.

Dataset | Duration | Child-speaking fraction
---|---|---
($mean\pm std$) |  ($mean\pm std$)
Pre-training | $14.05\pm 2.08$ | n.a
ADOSMod3 | $3.23\pm 1.61$ | $0.46\pm 0.18$
Simons | $19.05\pm 12.86$ | $0.40\pm 0.08$

##  4 System Description

###  4.1 Pre-training

Our research aims to adapt the existing self-supervised approaches to the
child-adult interaction domain through contrastive learning. Similar to [30],
our contrastive learning framework is based on the assumption that neighboring
segments from audio samples are highly likely to contain identical
information. For instance, it is probable that adjacent audio frames are
produced by the same speaker and are expected to contain similar semantic
meaning, linguistic content, as well as acoustic characteristics. To
elaborate, we define the dataset of audio samples as $N$, where each audio
sample is denoted as $x_{i}$. The corresponding neighboring audio segment is
represented as $x^{\prime}_{i}$, and is defined as any audio sample that has a
time shift of half a second or less from the original sample $x_{i}$.

As outlined in the previous section, transformer-based models first transform
the input speech sample $x$ to intermediate features ${z}$ using the feature
encoder $f(\cdot)$ on the basis of CNNs. Subsequently, the transformer encoder
$g(\cdot)$ maps the features ${z}$ to contextualized representations ${c}$.
Consequently, we can create similar pairs of contextualized representations
$c_{i}$ and $c^{\prime}_{i}$ from the neighboring audio segments $x_{i}$ and
$x^{\prime}_{i}$, with the remaining pairs being considered as negative pairs:

$\text{Positive Pairs}:c_{i}\approx c^{\prime}_{i}$ (1) $\text{Negative
Pairs}:c_{i}\neq c_{k}\text{ , }c_{i}\neq c^{\prime}_{k},\text{where }i\neq k$
(2)

Motivated by SimCLR [31], we apply the NTXent contrastive loss [32] as the
pretraining objective with the adult-child conversational corpora. Given the
temperature value $\tau$, the loss function $L_{NTXent}$ for the positive
audio pairs $x_{i}$ and $x^{\prime}_{i}$ within a batch of $B$ input audio is:

$\scriptsize-
log{\frac{exp(sim(c_{i},c^{\prime}_{i})/\tau)}{\sum^{B}_{\begin{subarray}{c}k=0\\\
k\neq
i\end{subarray}}exp(sim(z_{i},z_{k})/\tau)+\sum^{B}_{\begin{subarray}{c}k=0\end{subarray}}exp(sim(z_{i},z^{\prime}_{k})/\tau)}}$
(3)

###  4.2 Downstream Classifier Architectures

We use two different neural network models for child-adult speaker
classification based on [33]. Both the classifiers include a self-attention
based projector module, whereas one of them uses CNNs to capture speaker
characteristics and the other uses Recurrent Neural Netowrks (RNN) to model
the temporal dependencies present in the signal.

The RNN-based classifier consists of a stacked sequence of a Feed Forward
Layer (FFL), a bidirectional Long Short Term Memory (LSTM) layer, a self-
attention based projector layer and an output layer comprised of 2 FFLs,
separated by a non-linear activation. The CNN classifier architecture is
comprised of a weighted feature extraction module, followed by a convolutional
module having 3 1D convolutional layers, each with a dropout and a non-linear
activation in between, a self-attention based projector layer and an output
layer comprised of 2 FFLs, separated by a non-linear activation. For all the
experiments we use Rectified Linear Unit (ReLU) as the non-linear activation
and a dropout ratio of $0.3$.

##  5 Experimental Setup

###  5.1 Child Adult Classification

In this study, we hypothesize leveraging unlabelled child-speech for pre-
training can guide models to learn the heterogeneous child speech and
interaction patterns, leading to enhanced performance of downstream child-
adult speaker classification. Instead of training from scratch, we pre-train
the existing W2V2 and WavLM models with additional unlabelled child speech by
unfreezing and updating specific transformer layers using a contrastive loss
described in Sec 4.1. We report the child-adult classification macro F1-score
on two labeled child-adult interaction corpus described in section 3. We
report the results in Table 3 and Table 4, where the first row denotes
downstream child-adult classification performance using the model solely
relying on pre-trained embeddings. The subsequent rows denote downstream task
performance using the models pre-trained with additional child speech, where
the number indicates the number of trainable transformer layers involved in
the pre-training task.

Table 2: Number of trainable parameters for the pre-training experiments based
on unfrozen transformer layers

Number of unfrozen transformer layers
---
1 | 2 | 3 | 4 | 5
6.2M | 13.5M | 20.8M | 27.1M | 33.8M

###  5.2 ADOSMod3 Experiments on Demographics

Table 3: Child-adult classification F1 score using W2V2. PT corresponds to
pre-training and the following number represents the number of layers used for
pre-training.

Model | ADOSMod3 | Simons
---|---|---
RNN | CNN | RNN | CNN
W2V2 - Base | 67.92 | 70.59 | 63.41 | 64.13
W2V2 - PT1 | 69.31 | 72.41 | 65.19 | 66.28
W2V2 - PT2 | 71.55 | 72.95 | 65.87 | 65.12
W2V2 - PT3 | 72.23 | 74.38 | 68.81 | 65.44
W2V2 - PT4 | 74.01 | 74.89 | 67.63 | 66.79
W2V2 - PT5 | 72.19 | 74.05 | 65.01 | 65.39

Table 4: Child-adult classification F1 score using WavLM pre-training. PT
corresponds to pre-training and the following number represents the number of
layers used for pre-training.

Model | ADOSMod3 | Simons
---|---|---
RNN | CNN | RNN | CNN
WavLM-Base | 72.73 | 73.09 | 71.78 | 70.25
WavLM - PT1 | 74.29 | 74.93 | 72.64 | 71.11
WavLM - PT2 | 76.66 | 75.81 | 72.88 | 72.74
WavLM - PT3 | 75.95 | 76.37 | 72.31 | 71.09
WavLM - PT4 | 75.18 | 75.92 | 72.01 | 71.59
WavLM - PT5 | 75.48 | 73.17 | 71.47 | 70.17

(a) Session A: W2V2-Base

(b) Session A: W2V2-PT4

(c) Session B: WavLM-Base

(d) Session B: WavLM-PT2

Figure 2: t-SNE plots of the most discriminative 2 components of the embedding
space corresponding to the classes

Our study also investigates the model performance across age-groups in the
ADOSMod3 corpus. Prior works [11, 8] have reported age as an important
variability factor impacting speech characteristics. Based on this hypothesis,
we conduct an experiment by partitioning the ADOSMod3 corpus ($3-13yrs$) into
three different age-groups (Age-group 1: 43-90 months, Age-group 2: 91-118
months, Age-group 3: 119-158 months), such that each group contains equal
number of sessions. For each of these groups, we report the child-adult
classification F1 score using the pre-trained base models of W2V2 and WavLM
and also the best-performing pre-trained models of those two categories.

Not only age, analyses of developmental changes in speech have revealed sex
("gender") differences in speech characteristics, especially post puberty [8].
In this work, we also report gender-based child-adult speaker classification
performance on ADOSMod3 dataset, with recordings from $244$ male and $84$
female individuals. Similar to age-focused experiments, the dataset is
partitioned into male and female subsets, and comparisons are drawn between
the base model and the best-performing pre-trained models for both W2V2 and
WavLM.

###  5.3 Experimental details

For both the pre-training and fine-tuning experiments, $Adam$ optimizer is
used with a batch size of $32$ samples and temperature is set to $0.1$. The
number of tunable parameters for the pre-training experiments is reported in
Table 2. The initial learning rate is set to 1e-5 and the models are trained
for $30$ epochs with an early stopping callback on validation loss, patience
being $5$ epochs. For the downstream child-adult classification task, the
model is trained to minimize the binary cross-entropy loss for a maximum of
$50$ epochs, while the initial learning rate for this experiment is 2e-4 with
a weight decay of 1e-4. For both the datasets, we use 70% for training, 15%
for validation and 15% for testing. We use the model checkpoints from
HuggingFace [34]. We pre-train the models using a single NVIDIA GeForce GPU
1080 Ti and each experiment took less than two days.

##  6 Results and Discussion

Figure 3: Gender based Child-adult classification F1 scores.

###  6.1 Classification Evaluation

In this subsection, we analyze the experimental results reported in Table 3
and Table 4 to address the following questions:

Does pre-training with more child speech improve the classification? The
results reveal that pre-training with additional child speech improves the
child-adult classification F1 score over the base model. This underscores the
models' ability to account for the heterogeneity that is inherent in
children's speech. It can be observed that WavLM-based pre-trained models show
better performance compared to W2V2 across all the experiments. Both the
classifiers show comparable performance, with the RNN-based classifier
yielding the best score in the majority of the experiments. Among the
datasets, the experimental results reveal better F1 scores in ADOSMod3
compared to the Simons corpus. One possible reason might be related to the
session recording length difference between the datasets. The average duration
of sessions in Simons corpus is much higher than ADOSMod3, resulting in
greater potential variability and heterogeneity, which may have degraded the
F1 scores.

Does pre-training with more transformer layers improve the classification? It
is interesting to note that, while in W2V2-based pre-training, the
classification F1 keeps improving by tuning more transformer layers, in the
case of WavLM, the performance improvements reach the maximum with tuning
fewer transformer layers. One possible explanation is that the WavLM model is
trained with an objective function to capture speaker related information,
helping the model to achieve the optimum performance with lesser training.
However, in both the scenarios, the model performance starts to degrade by
adding more than four transformer layers. As these models are designed to
provide generalized speech representations, tuning larger portions of these
pre-trained models on a relatively smaller dataset might lead to the loss of
generalizability, causing the performance to decrease for the classification
task. However, our results provide compelling evidence that it is beneficial
to adapt the last few transformer layers for the adult-child classification.

Qualitative analysis We present t-SNE visualizations of pre-trained embeddings
for $2$ output classes from two sessions in Figure 2. We plot the embeddings
with and without additional pre-training. In both cases, it is evident from
the plots that our method increases the discriminative information between
them.

###  6.2 Result Evaluation based on Demographics

For the gender-focused experiments, the relative improvement in F1 scores are
$6.39\%$ and $3.14\%$ for the male and female subsets. Possibly due to both
inherent speech pattern differences and inherent data distribution biases (see
Section 5.2), the models yield higher F1 scores in the male population than
the female population. For the age-focused experiment, the relative
improvements of $9.42\%$, $8.23\%$, and $4.06\%$ are seen for the three age-
groups (youngest to oldest). The results imply that it is intrinsically
challenging to model children of AG1 and AG2 due to the developing vocal tract
behaviors among these ages. As a consequence, adding more children speech in
training data provides greater benefits to the model to capture more relevant
information, resulting in more improvement in the younger age-group than the
older ones.

Figure 4: Age based Child-adult classification F1 scores.

##  7 Conclusion

Past work has demonstrated the promise of deploying self-supervised algorithms
in a variety of downstream tasks like ASR, speaker diarization, and speaker
verification [35, 19]. In this work, we investigate the utility of additional
pre-training with more child speech, even in the presence of the inherent
heterogeneity and variability, to improve child-adult speaker classification
in clinical recordings involving interactions with children with autism. The
experimental results with the proposed models support our hypothesis of
benefiting from incorporating child speech based additional pre-training,
across both age and gender dimensions of variability.

In this work, we used the manually-annotated ground truth labels for
identifying and evaluating the speech and non-speech regions. In the future,
we plan to build a child-adult diarization framework with an integrated Voice
Activity Detection (VAD) system to further reduce the need of human effort. In
addition, we plan to extend this study with an additional emphasis on early
vocalization and speech (from toddlers and infants) in the interaction. Unlike
verbally fluent children, toddler speech contains significant amounts of pre-
verbal sounds and non-verbal vocalizations, which pose additional challenges
for automated processing.

##  8 Acknowledgements

This work is supported by funds from USC Hearing, Communication and
Neuroscience (HCN) pre-doctoral fellowship, NIH and Simons foundation.

## References

  * [1] J. Volden and C. Lord, ``Neologisms and idiosyncratic language in autistic speakers,'' _Journal of autism and developmental disorders_ , vol. 21, no. 2, pp. 109–130, 1991.
  * [2] S. V. Huemer and V. Mann, ``A comprehensive profile of decoding and comprehension in autism spectrum disorders,'' _Journal of Autism and Developmental Disorders_ , vol. 40, pp. 485–493, 2010.
  * [3] S. H. Kim, R. Paul, H. Tager-Flusberg, and C. Lord, ``Language and communication in autism,'' _Handbook of Autism and Pervasive Developmental Disorders, Fourth Edition_ , 2014.
  * [4] T. Sorensen, E. Zane, T. Feng, S. Narayanan, and R. Grossman, ``Cross-modal coordination of face-directed gaze and emotional speech production in school-aged children and adolescents with asd,'' _Scientific reports_ , vol. 9, no. 1, p. 18301, 2019.
  * [5] Centers for Disease Control and Prevention (CDC), ``Mental health in the united states: parental report of diagnosed autism in children aged 4-17 years–united states, 2003-2004,'' _MMWR. Morbidity and mortality weekly report_ , vol. 55, no. 17, pp. 481–486, 2006.
  * [6] C. Lord, S. Risi, L. Lambrecht, E. H. Cook, B. L. Leventhal, P. C. DiLavore, A. Pickles, and M. Rutter, ``The autism diagnostic observation schedule—generic: A standard measure of social and communication deficits associated with the spectrum of autism,'' _Journal of autism and developmental disorders_ , vol. 30, pp. 205–223, 2000.
  * [7] R. Grzadzinski, T. Carr, C. Colombi, K. McGuire, S. Dufek, A. Pickles, and C. Lord, ``Measuring changes in social communication behaviors: preliminary development of the brief observation of social communication change (boscc),'' _Journal of autism and developmental disorders_ , vol. 46, pp. 2464–2479, 2016.
  * [8] S. Lee, A. Potamianos, and S. Narayanan, ``Acoustics of children’s speech: Developmental changes of temporal and spectral parameters,'' _The Journal of the Acoustical Society of America_ , vol. 105, no. 3, pp. 1455–1468, 1999.
  * [9] V. Bhardwaj, M. T. Ben Othman, V. Kukreja, Y. Belkhier, M. Bajaj, B. S. Goud, A. U. Rehman, M. Shafiq, and H. Hamam, ``Automatic speech recognition (asr) systems for children: A systematic literature review,'' _Applied Sciences_ , vol. 12, no. 9, p. 4419, 2022.
  * [10] P. Gurunath Shivakumar and S. Narayanan, ``End-to-end neural systems for automatic children speech recognition: An empirical study,'' _Computer Speech & Language_, vol. 72, p. 101289, 2022.
  * [11] R. Lahiri, M. Kumar, S. Bishop, and S. Narayanan, ``Learning domain invariant representations for child-adult classification from speech,'' in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 6749–6753.
  * [12] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, ``Bert: Pre-training of deep bidirectional transformers for language understanding,'' _arXiv preprint arXiv:1810.04805_ , 2018.
  * [13] A. T. Liu, S.-W. Li, and H.-y. Lee, ``Tera: Self-supervised learning of transformer encoder representation for speech,'' _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 29, pp. 2351–2366, 2021.
  * [14] S. Ling and Y. Liu, ``Decoar 2.0: Deep contextualized acoustic representations with vector quantization,'' _arXiv preprint arXiv:2012.06659_ , 2020.
  * [15] A. Baevski, Y. Zhou, A. Mohamed, and M. Auli, ``wav2vec 2.0: A framework for self-supervised learning of speech representations,'' _Advances in neural information processing systems_ , vol. 33, pp. 12 449–12 460, 2020.
  * [16] S. Schneider, A. Baevski, R. Collobert, and M. Auli, ``wav2vec: Unsupervised pre-training for speech recognition,'' _arXiv preprint arXiv:1904.05862_ , 2019.
  * [17] S. Pascual, M. Ravanelli, J. Serra, A. Bonafonte, and Y. Bengio, ``Learning problem-agnostic speech representations from multiple self-supervised tasks,'' _arXiv preprint arXiv:1904.03416_ , 2019.
  * [18] M. Ravanelli, J. Zhong, S. Pascual, P. Swietojanski, J. Monteiro, J. Trmal, and Y. Bengio, ``Multi-task self-supervised learning for robust speech recognition,'' in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 6989–6993.
  * [19] Z. Fan, M. Li, S. Zhou, and B. Xu, ``Exploring wav2vec 2.0 on speaker verification and language identification,'' _arXiv preprint arXiv:2012.06185_ , 2020.
  * [20] S. Chen, C. Wang, Z. Chen, Y. Wu, S. Liu, Z. Chen, J. Li, N. Kanda, T. Yoshioka, X. Xiao _et al._ , ``Wavlm: Large-scale self-supervised pre-training for full stack speech processing,'' _IEEE Journal of Selected Topics in Signal Processing_ , vol. 16, no. 6, pp. 1505–1518, 2022.
  * [21] Y.-A. Chung, W.-N. Hsu, H. Tang, and J. Glass, ``An unsupervised autoregressive model for speech representation learning,'' _arXiv preprint arXiv:1904.03240_ , 2019.
  * [22] Y. Wang, A. Boumadane, and A. Heba, ``A fine-tuned wav2vec 2.0/hubert benchmark for speech emotion recognition, speaker verification and spoken language understanding,'' _arXiv preprint arXiv:2111.02735_ , 2021.
  * [23] D. Álvarez, S. Pascual, and A. Bonafonte, ``Problem-agnostic speech embeddings for multi-speaker text-to-speech with samplernn,'' _arXiv preprint arXiv:1906.00733_ , 2019.
  * [24] M. Najafian and J. H. Hansen, ``Speaker independent diarization for child language environment analysis using deep neural networks,'' in _2016 IEEE Spoken Language Technology Workshop (SLT)_. IEEE, 2016, pp. 114–120.
  * [25] A. Cristia, S. Ganesh, M. Casillas, and S. Ganapathy, ``Talker diarization in the wild: The case of child-centered daylong audio-recordings,'' in _Interspeech 2018_ , 2018, pp. 2583–2587.
  * [26] T. Zhou, W. Cai, X. Chen, X. Zou, S. Zhang, and M. Li, ``Speaker diarization system for autism children's real-life audio data,'' in _2016 10th International Symposium on Chinese Spoken Language Processing (ISCSLP)_. IEEE, 2016, pp. 1–5.
  * [27] J. Xie, L. P. Garcia-Perera, D. Povey, and S. Khudanpur, ``Multi-plda diarization on children's speech.'' in _Interspeech_ , 2019.
  * [28] N. R. Koluguri, M. Kumar, S. H. Kim, C. Lord, and S. Narayanan, ``Meta-learning for robust child-adult classification from speech,'' in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 8094–8098.
  * [29] S. Krishnamachari, M. Kumar, S. H. Kim, C. Lord, and S. Narayanan, ``Developing neural representations for robust child-adult diarization,'' in _2021 IEEE Spoken Language Technology Workshop (SLT)_. IEEE, 2021, pp. 590–597.
  * [30] V. Sachidananda, S.-Y. Tseng, E. Marchi, S. Kajarekar, and P. Georgiou, ``Calm: Contrastive aligned audio-language multirate and multimodal representations,'' _arXiv preprint arXiv:2202.03587_ , 2022.
  * [31] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, ``A simple framework for contrastive learning of visual representations,'' in _International conference on machine learning_. PMLR, 2020, pp. 1597–1607.
  * [32] K. Sohn, ``Improved deep metric learning with multi-class n-pair loss objective,'' _Advances in neural information processing systems_ , vol. 29, 2016.
  * [33] T. Feng, R. Hebbar, and S. Narayanan, ``Trustser: On the trustworthiness of fine-tuning pre-trained speech embeddings for speech emotion recognition,'' _arXiv preprint arXiv:2305.11229_ , 2023.
  * [34] T. Wolf, L. Debut, V. Sanh, J. Chaumond, C. Delangue, A. Moi, P. Cistac, T. Rault, R. Louf, M. Funtowicz _et al._ , ``Transformers: State-of-the-art natural language processing,'' in _Proceedings of the 2020 conference on empirical methods in natural language processing: system demonstrations_ , 2020, pp. 38–45.
  * [35] S.-w. Yang, P.-H. Chi, Y.-S. Chuang, C.-I. J. Lai, K. Lakhotia, Y. Y. Lin, A. T. Liu, J. Shi, X. Chang, G.-T. Lin _et al._ , ``Superb: Speech processing universal performance benchmark,'' _arXiv preprint arXiv:2105.01051_ , 2021.


# Numerical approximation of viscous contact problems applied to glacial
sliding

Gonzalo G. de Diego<EMAIL_ADDRESS>Patrick E. Farrell  Ian J.
Hewitt  Mathematical Institute, University of Oxford, Oxford, OX2 6GG, UK

###### Abstract

Viscous contact problems describe the time evolution of fluid flows in contact
with a surface from which they can detach. These problems are of particular
importance in glaciology, where they arise in the study of grounding lines and
subglacial cavities. In this work, we propose a novel numerical method for
solving viscous contact problems based on a mixed formulation with Lagrange
multipliers of a variational inequality involving the Stokes equation. The
advection equation for evolving the geometry of the domain occupied by the
fluid is then solved via a specially-built upwinding scheme, leading to a
robust and accurate algorithm for viscous contact problems. We then use this
numerical scheme to reconstruct friction laws for glacial sliding with
cavitation. We also study the effects of unsteady water pressures on sliding
with cavitation.

###### keywords:

##  1 Introduction

Viscous contact problems are time-dependent fluid flow problems in which the
fluid is in contact with a solid surface from which it can detach. Contact
problems of this type arise when modelling glacial ice flow, which is
typically treated as a viscous fluid flow (Schoof & Hewitt, 2013). On a large
scale, they are relevant to marine ice sheets with a grounding line (Schoof,
2007, 2011) and, on a smaller scale, to the formation of subglacial cavities
when the ice slides over bedrock undulations (Fowler, 1986; Schoof, 2005;
Gagliardini et al., 2007). These problems share a very similar mathematical
structure and are of great importance for understanding ice sheet dynamics and
predicting future sea level rise.

In this paper, we present a novel numerical approach for solving viscous
contact problems in an accurate and robust way. This method relies on the
formulation of the contact problem as a variational inequality. There exist
many different approaches to solving variational inequalities with finite
element methods. Of these, we find that solving the contact problem with a
piecewise constant Lagrange multiplier is particularly well suited since it
allows us to satisfy a discrete version of the contact boundary conditions
exactly in the discretisation. We find that this latter property of the
numerical scheme allows us to evolve the surface of the viscous flow in a
particularly robust way. To the best of our knowledge, viscous contact
problems have only previously been solved as variational inequalities in
Stubblefield et al. (2021), where a penalty method is used.

Although the numerical method presented in this paper is suitable for any
viscous contact problem, here we focus on subglacial cavitation and its
application to glacial sliding. Glacial sliding represents a fundamental
component of glacier dynamics and is one of the major uncertainties in ice
sheet modelling (Ritz et al., 2015). It involves a variety of physical
mechanisms among which is cavitation. Sliding with cavitation has been studied
theoretically by means of linearised approaches (Fowler, 1986; Schoof, 2005)
and these investigations are limited to relatively simple scenarios with two-
dimensional geometries, smooth bedrocks with small variations and steady
conditions. Numerical methods, on the other hand, allow the extension to more
complicated situations in a more straightforward manner (Gagliardini et al.,
2007; Helanow et al., 2020). In this work, we attempt to more fully exploit
the mathematical structure of the viscous contact problem, to improve the
accuracy and robustness of the numerical methods employed.

In Section 2 of this paper, we introduce our algorithm for solving viscous
contact problems, using the setup of the glacial cavitation problem to provide
a concrete context. We then present two applications of the algorithm. First,
in Section 3, we compute the steady sliding law for flow over a sinusoidal bed
for linear and nonlinear rheologies, expanding on the results in Gagliardini
et al. (2007). Furthermore, in this section we compare our results with those
obtained from the linearised theory to validate the algorithm. Then, in
Section 4, we explore the effects of unsteady water pressures on glacial
sliding by calculating the basal sliding velocities and cavity shapes under
oscillating water pressures.

##  2 Formulation and approximation of the problem of subglacial cavitation

We focus here on formulating the viscous contact problem of subglacial
cavitation and describing a finite element scheme to approximate it.
Subglacial cavitation is mathematically very similar to the problem of a
marine ice sheet with a grounding line. Therefore, the extension of the finite
element scheme presented here to grounding line problems requires minor
modifications, such as the inclusion of non-zero tangential stress (a friction
law). Note that the problem described here is the same, though with some
different notation, as the subglacial cavitation problem studied in
Gagliardini et al. (2007). A linearised version, assuming small amplitude of
the bed bump, is equivalent to the problem studied in Fowler (1986) and Schoof
(2005).

Subglacial cavitation occurs at the ice-bedrock interface, over length scales
corresponding to the size of the bedrock obstacles. These length scales are
generally several orders of magnitude smaller than those of the glacier. For
this reason, the computational domain $\Omega$ in which we model the formation
of cavities is a thin layer of ice of finite height located under a larger
mass of ice, see Figure 1. We assume the bedrock, and therefore also $\Omega$,
to be periodic in the horizontal direction. The upper boundary $\Gamma_{t}$
represents a fictional boundary separating $\Omega$ from the ice above. The
height of the bedrock is given by the function $b(x)$ and the height of the
lower boundary by $\theta(x,t)$. The lower boundary is divided into an
attached region $\Gamma_{a}$ where $\theta(x,t)=b(x)$, and a detached (i.e.
cavitated) region $\Gamma_{d}$ where $\theta(x,t)>b(x)$.

Figure 1: The periodic domain $\Omega$ on which we model the evolution of a
subglacial cavity.

The ice is generally modelled with the Stokes equations. If we denote the
velocity and pressure by $\boldsymbol{u}=(u,v)$ and $p$ respectively, these
can be written as

$\displaystyle-\nabla\cdot\left(2\eta\varepsilon(\boldsymbol{u})\right)+\nabla
p$ $\displaystyle=0,$ (1a) $\displaystyle\nabla\cdot\boldsymbol{u}$
$\displaystyle=0.$ (1b)

Here, $\eta=\eta(|\varepsilon(\boldsymbol{u})|)$ is the effective viscosity of
the ice, which is usually modelled with Glen’s law (Glen, 1958),

$\eta(|\varepsilon(\boldsymbol{u})|)=\frac{1}{2}A^{-1/n}\left(\frac{1}{2}|\varepsilon(\boldsymbol{u})|^{2}\right)^{\frac{1-n}{2n}},$
(2)

where $n\approx 3$ for ice and $A$ is a potentially temperature-dependent
parameter, but which we treat as constant. The tensor
$\varepsilon(\boldsymbol{u})$ represents the strain rate of the ice and is
given by the symmetric component of the velocity gradient,

$\varepsilon(\boldsymbol{u})=\frac{1}{2}(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{\top}).$
(3)

For a given velocity and pressure field, we define the stress tensor
$\sigma=\sigma(\boldsymbol{u},p)$ by

$\sigma=2\eta\varepsilon(\boldsymbol{u})-pI,$ (4)

where $I$ is the identity tensor field. If $\boldsymbol{n}$ denotes the
outwards-pointing normal vector to the boundary of $\Omega$, the normal and
tangential stresses at the boundary are defined by

$\sigma_{nn}=(\sigma\boldsymbol{n})\cdot\boldsymbol{n}\quad\text{and}\quad\boldsymbol{\sigma}_{nt}=\sigma\boldsymbol{n}-\sigma_{nn}\boldsymbol{n}.$
(5)

The Stokes equations (1) must be complemented with a set of boundary
conditions. Along the cavitated region $\Gamma_{d}$ of the lower boundary, we
assume the ice to be in contact with water at a pressure $p_{w}$, and we
prescribe

$\sigma_{nn}=-p_{w}\quad\text{and}\quad\boldsymbol{\sigma}_{nt}=0\quad\text{on
$\Gamma_{d}$}.$ (6)

On the attached region $\Gamma_{a}$, we assume the ice to be lubricated by a
thin layer of water connected to the subglacial drainage system. For this
reason, we allow the ice to slide freely and set $\boldsymbol{\sigma}_{nt}=0$.
The possibility of detachment is realised by enforcing the contact boundary
conditions

$\boldsymbol{u}\cdot\boldsymbol{n}\leq 0,\quad\sigma_{nn}\leq-
p_{w}\quad\text{and}\quad(\boldsymbol{u}\cdot\boldsymbol{n})(\sigma_{nn}+p_{w})=0\quad\text{on
$\Gamma_{a}$}.$ (7)

Notice that here we are enforcing the ice to either remain attached (i.e.
$\boldsymbol{u}\cdot\boldsymbol{n}=0$) whenever the compressive normal stress
is larger than the water pressure or to have the possibility of detaching
($\boldsymbol{u}\cdot\boldsymbol{n}<0$) if the stress equals the water
pressure.

On the top boundary $\Gamma_{t}$, we enforce the boundary condition
$\sigma_{nn}=-p_{i}$, where $p_{i}$ is the overburden ice pressure. Finally,
we close the system with either the Dirichlet boundary condition

$u=u_{i}\quad\text{on $\Gamma_{t}$},$ (8)

or the Neumann boundary condition

$\boldsymbol{\sigma}_{nt}=\tau_{b}\,(1,0)^{\top}\quad\text{on $\Gamma_{t}$},$
(9)

where $\tau_{b}$ is the basal shear stress. As explained further on in Section
3, we remark that the horizontal velocity $u_{i}$ is not the same as the basal
sliding speed (though they become approximately equal in the linearised
theory).

The Stokes equation must be coupled to a local advection equation which
describes the evolution in time of the cavity roof:

$\displaystyle\frac{\partial\theta}{\partial
t}+\boldsymbol{u}\cdot\left(\frac{\partial\theta}{\partial
x},-1\right)^{\top}$ $\displaystyle=0,$ (10a) $\displaystyle\theta$
$\displaystyle\geq b.$ (10b)

Since $\theta$ is time-dependent, so is the domain $\Omega$ and the attached
and detached regions $\Gamma_{a}$ and $\Gamma_{d}$ along the lower boundary.

Our goal in the rest of this section is to make use of the fact that the
instantaneous Stokes problem (1) together with the contact boundary conditions
(7) is equivalent to a variational inequality. This mathematical structure can
be exploited by a finite element discretisation.

###  2.1 A Stokes variational inequality

The finite element method can be used to compute approximate solutions of the
Stokes equation with the contact boundary conditions (7). As a first step, we
must find a weak formulation by multiplying (1) by suitable test functions and
integrating by parts. In order to keep the notation simpler, we assume in this
section and in Section 2.2 that we enforce the Neumann boundary condition (9)
instead of (8). In this case, a suitable space of test functions is the vector
Sobolev space $V=\boldsymbol{W}^{1,\kappa}(\Omega)$ for the velocity and the
Lebesgue space $Q=L^{\kappa^{\prime}}(\Omega)$ for the pressure. Here,
$\kappa=1+1/n$, where $n$ is the coefficient in Glen’s law, and
$\kappa^{\prime}=\kappa/(\kappa-1)$ is the Hölder conjugate of $\kappa$.

Due to the contact boundary conditions (7), the associated weak formulation of
the Stokes equation is a variational inequality. Given the convex set of
admissible velocity fields

$K=\left\\{\boldsymbol{v}\in V:\boldsymbol{v}\cdot\boldsymbol{n}\leq 0\text{
on $\Gamma_{a}$}\right\\},$ (11)

the weak formulation of (1) is: find $(\boldsymbol{u},p)\in K\times Q$ such
that

$\displaystyle
a(\boldsymbol{u},\boldsymbol{v}-\boldsymbol{u})-b(p,\boldsymbol{v}-\boldsymbol{u})$
$\displaystyle\geq f(\boldsymbol{v}-\boldsymbol{u})$
$\displaystyle\forall\boldsymbol{v}\in K,$ (12a) $\displaystyle
b(q,\boldsymbol{u})$ $\displaystyle=0$ $\displaystyle\forall q\in Q.$ (12b)

Here, the operators $a:V\times V\to\mathbb{R}$ and $b:Q\times V\to\mathbb{R}$
are defined by

$a(\boldsymbol{u},\boldsymbol{v})=\int_{\Omega}\eta\,\varepsilon(\boldsymbol{u}):\varepsilon(\boldsymbol{v})\,\mathrm{d}x\quad\text{and}\quad
b(q,\boldsymbol{v})=\int_{\Omega}q(\nabla\cdot\boldsymbol{v})\,\mathrm{d}x,$
(13)

and the linear functional $f:V\to\mathbb{R}$ by

$f(\boldsymbol{v})=\int_{\Gamma_{t}}\tau_{b}(\boldsymbol{v}\cdot(1,0))\,\mathrm{d}s-\int_{\Gamma_{t}}p_{i}(\boldsymbol{v}\cdot\boldsymbol{n})\,\mathrm{d}s-\int_{\Gamma_{a}\cup\Gamma_{d}}p_{w}(\boldsymbol{v}\cdot\boldsymbol{n})\,\mathrm{d}s.$
(14)

If the water pressure $p_{w}$ is constant along the bedrock, as we shall
assume, the incompressibility of the flow velocity can be used to rewrite the
right hand side as follows:

$f(\boldsymbol{v})=\int_{\Gamma_{t}}\tau_{b}(\boldsymbol{v}\cdot(1,0))\,\mathrm{d}s-\int_{\Gamma_{t}}N(\boldsymbol{v}\cdot\boldsymbol{n})\,\mathrm{d}s,$
(15)

where $N=p_{i}-p_{w}$ is the effective pressure. Consequently, the solution
only depends on $N$, not on the particular values of $p_{i}$ and $p_{w}$.

The variational inequality (12) is equivalent to the Stokes equation (1) in
the sense that if the velocity $\boldsymbol{u}$ is at least twice
differentiable and the pressure $p$ differentiable, then integration by parts
can be performed in order to arrive at (1) with the boundary conditions
introduced in the previous section. Of the three conditions in (7), only the
kinematic condition $\boldsymbol{u}\cdot\boldsymbol{n}\leq 0$ is enforced
explicitly in the definition of $K$; the remaining conditions are enforced
implicitly in (12).

There exist many different approaches for solving a variational inequality
with the finite element method (Glowinski et al., 1981). For example, in the
penalty method, the constraint $\boldsymbol{u}\cdot\boldsymbol{n}\leq 0$ in
the definition of $K$ is enforced via the addition of a nonlinear penalty term
in (12a) and the variational inequality is transformed into a variational
equality over $\boldsymbol{W}^{1,\kappa}(\Omega)$. This is the approach used
in Stubblefield et al. (2021). Of these approaches, we find that what works
best in combination with the advection equation (10) is a mixed method with
Lagrange multipliers.

The contact boundary conditions (7) can be explicitly enforced in (12) with a
Lagrange multiplier which essentially represents the sum of the normal stress
and the water pressure $\sigma_{nn}+p_{w}$ along the attached region
$\Gamma_{a}$. In order to introduce this reformulation of (12), we need to
define a suitable set of admissible Lagrange multipliers. We define the range
of the normal trace operator onto $\Gamma_{a}$ by

$\Sigma=\left\\{\boldsymbol{v}\cdot\boldsymbol{n}|_{\Gamma_{a}}:\boldsymbol{v}\in\boldsymbol{W}^{1,\kappa}(\Omega)\right\\}.$
(16)

The space of admissible Lagrange multipliers is then the convex cone of
elements in the dual space $\Sigma^{\prime}$ satisfying a weak equivalent of
the dynamic contact boundary condition $\sigma_{nn}+p_{w}\leq 0$ in (7):

$\Lambda=\left\\{\mu\in\Sigma^{\prime}:\mu(\xi)\geq
0\quad\text{$\forall\xi\in\Sigma$ such that $\xi\leq 0$ on
$\Gamma_{a}$}\right\\}.$ (17)

The mixed formulation of the variational inequality (12) may be written as:
find $(\boldsymbol{u},p,\lambda)\in V\times Q\times\Sigma^{\prime}$ such that

$\displaystyle
a(\boldsymbol{u},\boldsymbol{v})-b(p,\boldsymbol{v})-\lambda(\boldsymbol{v}\cdot\boldsymbol{n}|_{\Gamma_{a}})=f(\boldsymbol{v})$
$\displaystyle\forall\boldsymbol{v}\in V,$ (18a) $\displaystyle
b(q,\boldsymbol{u})=0$ $\displaystyle\forall q\in Q,$ (18b)
$\displaystyle\mu(\boldsymbol{u}\cdot\boldsymbol{n}|_{\Gamma_{a}})\geq
0\quad\forall\mu\in\Lambda,\quad\lambda\in\Lambda\quad\text{and}\quad\lambda(\boldsymbol{u}\cdot\boldsymbol{n}|_{\Gamma_{a}})=0$
. (18c)

Equation (18a) can be obtained by multiplying (1a) by a test function
$\boldsymbol{v}$, integrating by parts, and setting

$\lambda=\sigma_{nn}+p_{w}.$ (19)

Formulation (18) offers several advantages for finite element approximation in
comparison to (12). First, there is no longer any need to build an
approximation of the convex set $K$ because we seek the velocity in the space
$V$. Second, the conditions (18c), which are a weak analogue of the contact
boundary conditions (7), can be rewritten as a nonlinear equation at the
discrete level, see (22c) below. This means that the contact boundary
conditions can be enforced exactly by solving the finite element system
analogous to (18). As we explain shortly, this leads to a robust and stable
numerical method for solving the complete subglacial cavitation system.

Under certain conditions on the geometry of the domain and the linear
functional $f$, the variational inequality (12) has a unique solution. A proof
can be found in de Diego et al. (2021), together with an analysis of the
finite element approximation of the mixed system (18).

###  2.2 A numerical method for subglacial cavitation

In this section we present a numerical algorithm for solving the complete
subglacial cavitation problem which results from coupling the Stokes equation
(1) to the time-dependent advection equation (10). We write the discrete
counterparts to these equations and explain how they are coupled by deforming
the domain according to a contact criterion. The resulting algorithm is
summarised in Algorithm 1. The solver is implemented in Firedrake (Rathgeber
et al., 2016), using the version available at zenodo/Firedrake-20211103.0
(2021). The code for the viscous contact solver presented in Algorithm 1 is
openly available in the archived repository de Diego (2021) and in
https://bitbucket.org/gonzalogddiego/subglacialcavitysolver_2021/src/master/.

Subglacial cavitation is a time-dependent problem, so its discretisation
requires a partition of a given time interval $[0,T]$ into intervals of length
$\Delta t$. At each time step $t_{k}=k\Delta t$ for $k=0,1,...,N=T/\Delta t$,
we must consider a discrete cavity roof $\theta_{h}^{k}$ and domain
$\Omega^{k}$ given by

$\Omega^{k}=\left\\{(x,y)\in\mathbb{R}^{2}:0\leq x\leq
L,\quad\theta_{h}^{k}(x)\leq y\leq H\right\\}.$ (20)

Additionally, the lower boundary is the union of the attached and detached
regions $\Gamma_{a}^{k}$ and $\Gamma_{d}^{k}$ respectively. These are also
time-dependent and are determined at each time step by the contact criterion
(28) below.

Set $\theta^{0}$ ;

for  _$k=0,1,2,\dots$_ do

Set $\Gamma^{k}_{a}$ and $\Gamma^{k}_{d}$ according to (28);

Find $(\boldsymbol{u}^{k}_{h},p^{k}_{h},\lambda^{k}_{h})$ such that (22) holds
for $\theta=\theta^{k}$ and $\Gamma_{a}=\Gamma^{k}_{a}$;

Calculate $\theta^{k+1}$ with (27) ;

If $\theta_{i}^{k+1}<b(x_{i})$, set $\theta_{i}^{k+1}=b(x_{i})$;

Deform mesh;

end for

Algorithm 1 Solution procedure for the subglacial cavity problem

####  2.2.1 Finite element approximation of the Stokes variational inequality

On each domain $\Omega^{k}$, we seek a finite element approximation to (18).
To this end, we consider a triangulation of $\Omega^{k}$ and the finite
element spaces $V_{h}$, $Q_{h}$ and $\Sigma_{h}$ given by piecewise continuous
quadratic vector fields, piecewise constant scalar fields and piecewise
constant scalar functions on $\Gamma_{a}^{k}$ respectively. This velocity-
pressure pair is known to be inf-sup stable in two dimensions (Boffi et al.,
2013). Furthermore, as we explain in Section 2.2.2, the use of piecewise
constant elements for the Lagrange multiplier enables the use of a very simple
upwinding scheme for solving the advection equation (10) in a manner
consistent with the discrete contact boundary conditions.

We write $V_{h}=\operatorname*{span}{\\{\boldsymbol{v}_{i}\\}}_{i=1}^{N_{v}}$,
$Q_{h}=\operatorname*{span}{\\{q_{j}\\}}_{j=1}^{N_{q}}$ and
$\Sigma_{h}=\operatorname*{span}{\\{\mu_{k}\\}}_{k=1}^{N_{\mu}}$, where
$N_{v}=\dim{V_{h}}$, $N_{q}=\dim{Q_{h}}$ and $N_{\mu}=\dim{\Sigma_{h}}$. For
the functions $(\boldsymbol{u}_{h},p_{h},\lambda_{h})\in V_{h}\times
Q_{h}\times\Sigma_{h}$, we denote the respective degrees of freedom (DoFs) in
$\mathbb{R}^{N_{v}}$, $\mathbb{R}^{N_{q}}$ and $\mathbb{R}^{N_{\mu}}$ by
$\mathbf{u}$, $\mathbf{p}$ and $\boldsymbol{\lambda}$. In order to write an
algebraic counterpart of (18c), we need the operator

$\boldsymbol{\gamma_{n}}:\mathbb{R}^{N_{v}}\to\mathbb{R}^{N_{\mu}}$ (21)

that returns the average normal components of a vector $\boldsymbol{v}_{h}\in
V_{h}$ along the edges on $\Gamma_{a}^{k}$. The discrete counterpart to (18)
is

$\displaystyle\mathbf{A}(\mathbf{u})-\mathbf{B}\mathbf{p}-\mathbf{D}\boldsymbol{\lambda}$
$\displaystyle=\mathbf{f},$ (22a) $\displaystyle-\mathbf{B}^{\top}\mathbf{u}$
$\displaystyle=0,$ (22b)
$\displaystyle\boldsymbol{\lambda}+\mathbf{C}(\boldsymbol{\lambda},\mathbf{u})$
$\displaystyle=0.$ (22c)

Here, we have introduced the matrices $\mathbf{B}\in\mathbb{R}^{N_{v}\times
N_{q}}$ and $\mathbf{D}\in\mathbb{R}^{N_{v}\times N_{\mu}}$, the vector
$\mathbf{f}\in\mathbb{R}^{N_{v}}$ and the nonlinear operators
$\mathbf{A}:\mathbb{R}^{N_{v}}\to\mathbb{R}^{N_{v}}$ and
$\mathbf{C}:\mathbb{R}^{N_{\mu}}\times\mathbb{R}^{N_{v}}\to\mathbb{R}^{N_{\mu}}$.
The matrices are given by the elements
$\mathbf{B}_{ij}=b(\boldsymbol{v}_{i},q_{j})$ and
$\mathbf{D}_{ij}=\int_{\Gamma^{k}_{a}}\mu_{j}(\boldsymbol{v}_{i}\cdot\boldsymbol{n})\,\mathrm{d}s$
and the vector $\mathbf{f}$ by

$\mathbf{f}_{i}=-\int_{\Gamma_{t}}p_{i}(\boldsymbol{v}_{i}\cdot\boldsymbol{n})\mathrm{d}s-\int_{\Gamma_{a}\cup\Gamma_{d}}p_{w}(\boldsymbol{v}_{i}\cdot\boldsymbol{n})\mathrm{d}s+\int_{\Gamma_{t}}\tau_{b}(\boldsymbol{v}_{i}\cdot(1,0))\,\mathrm{d}s.$
(23)

The nonlinear operators are given by

$[\mathbf{A}(\mathbf{u})]_{i}=a(\boldsymbol{u}_{h},\boldsymbol{v}_{i})$ (24)

and

$\mathbf{C}(\boldsymbol{\lambda},\mathbf{u})=\max{\left\\{0,-\boldsymbol{\lambda}+c(\boldsymbol{\gamma_{n}}\mathbf{u})\right\\}}$
(25)

for an arbitrary $c>0$. Here, the $\max$ operation is understood to be carried
out for each of the elements in the vector
$-\boldsymbol{\lambda}+c(\boldsymbol{\gamma_{n}}\mathbf{u})\in\mathbb{R}^{N_{\mu}}$.
The use of (25) is a common way of expressing contact conditions; a particular
advantage is that the nonlinear system (22) can be solved with a semi-smooth
Newton method (Hintermüller et al., 2002). Moreover, (22c) is equivalent to
enforcing

$\boldsymbol{\gamma_{n}}\mathbf{u}\leq 0,\quad\boldsymbol{\lambda}\leq
0\quad\text{and}\quad(\boldsymbol{\gamma_{n}}\mathbf{u})\cdot\boldsymbol{\lambda}=0,$
(26)

a discrete equivalent of (7). By solving (26), we see that on each edge on
$\Gamma_{a}^{k}$ we either have that the edge wishes to detach
($\boldsymbol{\gamma_{n}}\mathbf{u}<0$) or to remain attached to the bed
($\boldsymbol{\lambda}\leq 0$ and $\boldsymbol{\gamma_{n}}\mathbf{u}=0$).

####  2.2.2 The discrete advection equation and a contact criterion

Figure 2: Notation and ordering for nodes $\\{x_{i}\\}$ and edges
$\\{e_{i}\\}$ along the attached region $\Gamma_{a}^{k}$.

In order to write a discretisation of (10), we introduce some further
notation. We denote the points along the lower boundary by
$(x_{i},\theta_{i})$ with the index $i$ increasing from left to right as in
Figure 2. The edge between $(x_{i-1},\theta_{i-1})$ and $(x_{i},\theta_{i})$
is denoted by $e_{i}$. We write $u^{k}_{n,i}$ for the value of
$\boldsymbol{\gamma_{n}}\mathbf{u}$ associated to $e_{i}$ at time $t_{k}$ (by
the definition of $\boldsymbol{\gamma_{n}}$, this value corresponds to the
average value of $\boldsymbol{u}\cdot\boldsymbol{n}$ along $e_{i}$).

A numerical scheme that solves the advection equation (10a) should satisfy two
conditions. First, it must be stable, as described in standard textbooks
(LeVeque, 2007). Additionally, the scheme should be compatible with the
discrete contact boundary conditions (26), in the sense that if for an edge
$e_{i}$ we have that $u_{n,i}^{k}=0$, then we should have that both
$\theta_{i-1}$ and $\theta_{i}$ remain unchanged. A naive discretisation of
(10) will detach all edges within a few timesteps due to approximation errors.
Here, we propose the following explicit scheme:

$\frac{\theta_{i}^{k+1}-\theta_{i}^{k}}{\Delta
t}=\left(\left(\dfrac{\theta^{k}_{i}-\theta^{k}_{i-1}}{x_{i}-x_{i-1}}\right)^{2}+1\right)^{1/2}u^{k}_{n,i},$
(27)

where the velocity at each node is taken from the edge located immediately
upstream. This results in an upwinding scheme with the property that
penetration cannot occur along $\Gamma_{a}$, because the average values of
$\boldsymbol{u}_{h}\cdot\boldsymbol{n}$ along the edges are non-positive, see
(26). However, a node on a detached edge can penetrate the bedrock. To avoid
this and enforce (10b), we simply set $\theta^{k}_{i}=b(x_{i})$ whenever
$\theta^{k}_{i}<b(x_{i})$ .

Finally, given a cavity roof $\theta_{h}^{k}$ at $t_{k}$, we must set the
attached and detached regions $\Gamma_{a}^{k}$ and $\Gamma_{d}^{k}$ to solve
the Stokes variational inequality at $t_{k}$. This is carried out by assigning
each edge $e_{i}$ with endpoints $(x_{i-1},\theta^{k}_{i-1})$ and
$(x_{i},\theta^{k}_{i})$ to either $\Gamma_{a}^{k}$ or $\Gamma_{d}^{k}$. We do
this by looking at the node downstream of the edge, in accordance with the
scheme for the advection equation:

$e_{i}$ belongs to $\Gamma_{d}^{k}$ if and only if $\theta^{k}_{i}>b(x_{i})$.
(28)

Consequently, given a cavity roof $\theta^{k}_{h}$ at time $t_{k}$, the
attached and detached regions $\Gamma^{k}_{a}$ and $\Gamma^{k}_{d}$ are fully
determined.

##  3 Steady sliding with cavitation

The sliding of a glacier over its bedrock has been widely studied since
Weertman’s seminal work in 1957 (Weertman, 1957). In general, these studies
attempt to build a function known as the sliding law that captures the steady
relationship between the basal sliding speed $u_{b}$, the basal shear stress
$\tau_{b}$ and other variables such as the water pressure $p_{w}$. This
sliding law can then be used to prescribe a boundary condition at the ice-
bedrock interface in large-scale glacier models which do not resolve the
smaller-scale shape of that interface.

As a first application of the algorithm described in the previous section, we
build the sliding law for ice flowing over a sinusoidal bed. We first present
detailed results for a particular steady state in Section 3.1 to evaluate the
accuracy of the solver and the effect of mesh refinement. Then, in Section 3.2
we compute sliding laws for different values of the parameter $n$ in Glen’s
law (2).

The computations presented in this section and in Section 4 have been carried
out under the assumption that the water pressure $p_{w}$ is uniform along the
lower boundary. As explained in the previous section, the problem then only
depends on the effective pressure $N=p_{i}-p_{w}$, rather than individually on
$p_{i}$ or $p_{w}$ (the boundary conditions become $\sigma_{nn}=-N$ on
$\Gamma_{t}$, $\sigma_{nn}=0$ on $\Gamma_{d}$ and $\sigma_{nn}\leq 0$ on
$\Gamma_{a}$). We non-dimensionalise the calculations by scaling lengths with
the wavelength of the bed $L$, velocities with a characteristic velocity scale
$U$, and stresses with the characteristic stress scale $(U/(2AL))^{1/n}$. The
resulting problem depends only on the non-dimensional amplitude $r$ of the
sinusoidal bed, the non-dimensional velocity $u_{i}$ prescribed along the
upper boundary, and the non-dimensional effective pressure $N$.

###  3.1 Steady subglacial cavities

Table 1: Information about the different meshes used to compute the steady cavity states together with calculations of the basal shear stress $\tau_{b}$, the basal sliding speed $u_{b}$, and the detachment and reattachment points of the cavity. $n_{e}$ | mesh cells | $\tau_{b}$ | $u_{b}$ | detachment | reattachment
---|---|---|---|---|---
16 | 96 | 0.014772 | 0.98667 | 0.7500 | 1.0000
32 | 192 | 0.015143 | 0.98633 | 0.7188 | 1.0000
64 | 768 | 0.015484 | 0.98598 | 0.7188 | 1.0000
128 | 3072 | 0.015679 | 0.98577 | 0.7109 | 1.0000
192 | 7296 | 0.015741 | 0.98570 | 0.7135 | 0.9948
Figure 3: The steady cavity shape and normal stresses along the attached
region computed with different mesh sizes. In black, the stress distribution
obtained from the linearised theory. Here, the amplitude of the obstacle is
set to $r=0.01$, the scaled effective pressure to $N=0.3$, and the scaled
horizontal velocity at the top to $u_{i}=1$.

The steady state of equations (1) and (10) with the boundary conditions
$u=u_{i}$ and $\sigma_{nn}=-N$ on the top boundary $\Gamma_{t}$ can be found
with Algorithm 1 by evolving the cavity from an initial state until the norm
of the discrete derivative $(\theta^{k+1}-\theta^{k})/\Delta t$ is below a
prescribed threshold of $10^{-4}$. In this section, we find the steady cavity
for a Newtonian flow ($n=1$ in Glen’s law (2)) over a sinusoidal bed of small
amplitude $r=0.01$ given by $b(x)=rL\sin{\left(\pi/2+2\pi x/L\right)}$. The
non-dimensional effective pressure is set to $N=0.3$ and the non-dimensional
horizontal velocity to $u_{i}=1$. We consider a linear rheology and a bedrock
of small amplitude in order to compare our results with the analytical
solution of the linearised cavitation problem considered in Fowler (1986) and
Schoof (2005).

We use five different meshes with $n_{e}$ vertices uniformly distributed along
the lower boundary. In Table 1 we present the non-dimensional basal shear
stress $\tau_{b}$ and sliding speed $u_{b}$ along the cavities. These values
are calculated with

$\tau_{b}=-\frac{1}{L}\int_{\Gamma_{a}\cup\Gamma_{d}}(\sigma_{nn}+p_{w})\,n_{x}\,\mathrm{d}s\quad\text{and}\quad
u_{b}=\frac{1}{L}\int_{\Gamma_{a}\cup\Gamma_{d}}u\,\mathrm{d}s,$ (29)

where $\boldsymbol{n}=(n_{x},n_{y})$. Recall that the Lagrange multiplier
$\lambda$ represents the normal stress $\sigma_{nn}+p_{w}$ along $\Gamma_{a}$,
see (19); therefore $\tau_{b}$ can be calculated from $\lambda$ via

$\tau_{b}=-\frac{1}{L}\int_{\Gamma_{a}}\lambda n_{x}\,\mathrm{d}s.$ (30)

As we show in the next section, the values of $\tau_{b}$ and $u_{b}$ can be
used to construct a sliding law.

The formula for $u_{b}$ presented in (29) might seem strange if the subglacial
cavity domain is interpreted as a boundary layer between an ice sheet and the
bedrock. This would suggest we take the sliding speed to be the average value
of the horizontal velocity along the top boundary $\Gamma_{t}$. However, our
computations indicate that, if the height of the domain $H$ is sufficiently
large, we can expect the shear stress to approach a constant value and the
horizontal velocity to vary with $y^{n}$ as $y$ approaches $H$, where $n$ here
is the exponent in Glen’s law (2). Therefore, the horizontal velocity along
the top boundary depends strongly on the height of the domain. For this
reason, and following Gagliardini et al. (2007), we use the equation in (29)
to calculate $u_{b}$. In this case, we find that $u_{b}$ is independent of $H$
for sufficiently large values of $H$. In particular, throughout this paper we
set $H=L$. In agreement with Gagliardini et al. (2007), we find this value of
$H$ to be sufficiently large.

In Figure 3 we present the steady cavity shape and normal stresses
$\sigma_{nn}$ along the attached region $\Gamma_{a}$ for the different meshes.
We can see from these figures that the cavity shape is accurately computed
even with the coarsest mesh. Additionally, we also present the stress
distribution obtained from the linearised theory, which is uniquely determined
for an effective pressure $N$ and a sliding speed $u_{b}$. The result from the
linearised theory plotted in Figure 3 is computed with the value of $u_{b}$
calculated with the most refined mesh.

The plot for the normal stress distribution demonstrates that the contact
conditions (7) are satisfied exactly at the discrete level for all of the
meshes, because $\sigma_{nn}+p_{w}\leq 0$. This plot also exhibits the
singularity of the normal stresses at the reattachment point. This singularity
complicates the approximation of the normal stresses along the attached region
and can lead to very inaccurate computations of the sliding law in largely
cavitated states. However, Figure 3 also indicates that, with increasing mesh
refinement, the solver appears to converge towards the linearised solution.

At the discrete level, the Lagrange multiplier $\lambda$, which we use to
calculate $\tau_{b}$, is piecewise constant on each edge along the lower
boundary. In Figure 3, these values, the DoFs of $\lambda$, are plotted at the
midpoints of each edge. Due to the contact criterion (28), the edge
immediately upstream of the first reattached node is treated as part of
$\Gamma_{a}$ (this is required in order to allow for the possibility of
subsequent detachment there), which explains why there are non-zero values of
$\sigma_{nn}+p_{w}$ left of the reattachment node in Figure 3. This is
particularly visible for the coarsest mesh with $n_{e}=16$.

###  3.2 Computation of the linear and nonlinear steady sliding law

We next perform similar calculations to those of Section 3.1 but for varying
effective pressure $N$ and power-law exponent $n$. This allows us to map out a
steady sliding law for ice sliding over a hard bed with cavitation as in
Gagliardini et al. (2007). The steady sliding law can be mapped out by only
varying $N$ and keeping $u_{i}$ fixed, since dimensional analysis of the
steady problem shows that the scaled basal shear stress $\tau_{b}/N$ depends
only on the ratio $u_{b}/(ALN^{n})$, and not independently on $u_{b}$ or $N$
(Fowler, 1986) (the same will not be true of the unsteady problem in Section
4).

Several previous studies have suggested what form the sliding law should take,
both with and without cavitation (Kamb, 1970; Fowler, 1981, 1986; Gudmundsson,
1997; Schoof, 2005; Gagliardini et al., 2007). The law proposed in Gudmundsson
(1997) for the uncavitated case can be written as

$\left(\frac{\tau_{b}}{rN}\right)^{n}=\alpha(n)\frac{r}{AL}\frac{u_{b}}{N^{n}},$
(31)

where $r=a/L$ and $\alpha(n)$ is a function depending on $n$. The function
$\alpha(n)$ is related to the parameter $c_{0}$ (which also depends on $n$)
considered in Gudmundsson (1997) by

$\alpha(n)=\frac{(2\pi)^{n+2}}{2c_{0}}.$ (32)

For a Newtonian flow, the complex analysis method presented in Fowler (1986)
and Schoof (2005) yields an exact solution to the linearised problem. In
particular, for high effective pressures, no cavitation occurs and a linear
sliding law as in (31) with $c_{0}=1$ is found. For effective pressures larger
than a critical value $8\pi^{2}r\eta u_{b}/L$, cavitation occurs and the
sliding law becomes non-linear, varying with $N$ as well as $u_{b}$.

Figure 4: Computed sliding law for steady glacial sliding with cavitation
(left) and cavity endpoints (right). The parameter $n$ is Glen’s law (2) is
set to $n=1$ (green), 3 (blue) and 5 (red); for each $n$, we compute the
sliding law for beds of amplitude $r=0.01$ (circles, light), 0.04 (squares)
and 0.08 (triangles, dark). For these computations, we set the scaled velocity
at the top of the domain to $u_{i}=1$ and the scaled viscosity parameter to
$A=0.5$. We use a mesh with 192 cells along the lower boundary. The parameter
$\alpha(n)$ is computed from the slope of the curve near the origin for the
lowest $r$.

We compute the sliding law over a sinusoidal bed of different amplitudes $r$
and for $n=1$, 3 and 5, and we plot the results in Figure 4 using the scaling
suggested by (31). The mesh has 192 vertices along the lower boundary and 7296
cells. The location of the cavity endpoints is also plotted in Figure 4, along
with the solution to the linearised problem calculated with the method from
Fowler (1986) and Schoof (2005). For each $n$, the parameter $\alpha(n)$ is
computed by calculating the slope of the curve near the origin (where there is
no cavitation) for the lowest value of $r$. The corresponding values of the
parameter $c_{0}$ can then be calculated from (32); these values can be found
in Table 2 together with those obtained in Gudmundsson (1997) and Gagliardini
et al. (2007).

The computed sliding laws with cavitation in Figure 4 are multivalued for
$\tau_{b}/(rN)$ as expected (Fowler, 1986; Schoof, 2005). This aspect of the
law justifies the use of the Dirichlet boundary condition $u=u_{i}$ instead of
the alternative Neumann boundary condition. We find that if we use the Neumann
boundary condition $\boldsymbol{\sigma}_{nt}=\tau_{b}(1,0)$ and initiate the
cavity from a fully attached state, the solver always evolves to the steady
state associated to the upsloping region of the curve (see also Section 4).
From Figure 4, we also deduce that the validity of the sliding law (31) along
the linear segment of the curves (where either little or no cavitation has
occurred) decreases with increasing values of $r$ and $n$. For example, when
$n=5$, one can observe that the linear segment of the sliding law for $r=0.08$
clearly does not collapse onto the corresponding linear segment for $r=0.01$.
As soon as the cavity size increases and the sliding laws cease to be linear,
the aspect of these curves largely differ for different values of $n$.

In Figure 4, we use a different scaling to the one used in Gagliardini et al.
(2007). In Gagliardini et al. (2007), the computed maximum value reached by
$\tau_{b}/N$ is included in the scaling for the sliding law. In this way, the
maximum value reached by the scaled sliding law equals 1 by design. However,
we preferred the scaling based on (31) because it contains fewer terms that
are unknown a priori. It is also worth mentioning that, for different values
of $n$, the curves in Figure 4 do not collapse into a single curve when
plotted with the scaling from Gagliardini et al. (2007).

For the linear case with $n=1$, the numerical results computed with the finite
element solver highly resemble those obtained with the linearised solution.
For $r=0.08$ a slight difference with the linearised solution can be seen near
the peak of the sliding law. This difference is probably a consequence of
nonlinear effects that are accentuated with increasing amplitudes of bedrock
roughness.

Table 2: Value of the parameter $c_{0}$ associated to the sliding laws. This parameter is associated to the slope of these curves near the origin. | $r=0.01$ | Gudmundsson (1997) | Gagliardini et al. (2007)
---|---|---|---
$n=1$ | 1.0014 | 0.9936 | 0.9771
$n=3$ | 0.3434 | 0.3294 | 0.2769
$n=5$ | 0.1255 | 0.1153 | -

##  4 Unsteady sliding

In the previous section, the sliding law was constructed by computing steady
cavity states. However, field measurements from alpine glaciers and from the
Greenland Ice Sheet have found short term variations in the water pressure, on
timescales down to hours (Iken, 1981; Iken & Bindschadler, 1986; Sugiyama &
Gudmundsson, 2004; Andrews et al., 2014; Hoffman et al., 2016). In these
studies, variations of water pressure have been correlated with variations in
surface speeds, vertical strain and uplift. Subglacial cavitation has been
considered a possible mechanism causing these correlations (Iken &
Bindschadler, 1986; Mair et al., 2002; Sugiyama & Gudmundsson, 2004). These
observations motivate an investigation of glacier sliding under unsteady
conditions. In this section, we therefore compute the evolution in time of
subglacial cavities under oscillating water pressures and calculate the
corresponding unsteady basal sliding speeds and shear stresses. The study
published by Iken of the transient stages between steady cavity shapes (Iken,
1981) is the only numerical investigation of unsteady cavitation known to the
authors.

We initialise the computations from a steady state corresponding to a point in
the sliding law determined by an effective pressure $N_{0}$, a basal sliding
speed $u_{b,0}$, and a basal shear stress $\tau_{b,0}$. Instead of prescribing
the Dirichlet boundary condition $u=u_{i}$ on $\Gamma_{t}$, we enforce the
Neumann boundary condition $\boldsymbol{\sigma}_{nt}=\tau_{b,0}(1,0)$ on
$\Gamma_{t}$. We consider it more physically realistic to have the basal shear
stress fixed rather than the sliding speed because we can expect the basal
stresses to balance the gravitational driving stresses, which are essentially
fixed on these timescales (in practice, if water pressure variations are
spatially localised, the driving stress can be transferred to neighbouring
regions of the bed, but it is not easy to account for this within the current
boundary-layer treatment of the problem). We set $n=3$ in Glen’s law to model
the nonlinear rheology of ice. Algorithm 1 is used with a mesh with 192
elements along the lower boundary over a sinusoidal bed of amplitude $r=0.08$.

Figure 5: Evolution of the cavity volume $V$ (right axis) and basal sliding
speed $u_{b}$ (left axis, also for $N$) for imposed oscillating effective
pressures $N$ and fixed basal shear stress $\tau_{b}$ around the states with
$u_{b,0}/(2ALN_{0}^{n})=0.0448$ (top), $0.0761$ (middle) and $0.6074$ (lower).
One period of each solution is superimposed on the steady sliding law in the
right hand panels, as indicated by the coloured dots

The effects of unsteady water pressures differ depending on the initial steady
state from which we evolve the cavity. To illustrate this, three different
points along the steady sliding law are evolved by oscillating the effective
pressure with an amplitude of $0.1N_{0}$ and a fixed non-dimensional frequency
of 0.4 (note that one non-dimensional time unit is approximately the time
taken for ice at the top of the domain to traverse one wavelength of the bed).
The results are plotted in Figure 5. In the top and middle panels, the steady
states are points on the upsloping section of the curve. These results
indicate that, with increasing cavitation, the amplitude of the sliding speed
increases. Moreover, the phase difference between the oscillating sliding
speeds and the cavity volumes also changes: in the top figure, one can observe
that the maximum sliding speed is reached when the cavity is still growing,
while in middle one it is reached at the time of maximum cavitation. On the
other hand, maximum sliding speeds occur at the times of lowest effective
pressure in all cases. Field measurements have also found maximum surface
speeds to take place at moments of maximum water pressures (Iken &
Bindschadler, 1986; Sugiyama & Gudmundsson, 2004). There are slight
oscillations in the computed sliding speed when the cavity volume is at its
largest. These are numerical artefacts due to the stress singularity at the
reattachment point of the cavity having an increasing effect on the overall
solution of the problem as the cavity volume grows. In these situations, a
small displacement of the reattachment point has a large effect on the stress
distribution along the bed and therefore also on the computed sliding speed.

The downsloping section of the curve produces a so-called rate-weakening
sliding regime in which, for a supposed fixed effective pressure, an increase
in the sliding speed is accompanied by a decrease in the basal drag. Rate-
weakening sliding has been observed in a laboratory setting for ice sliding
over a sinusoidal bed (Zoet & Iverson, 2015), although it has been questioned
whether such a sliding regime can arise for more realistic bed geometries
(Fowler, 1987; Schoof, 2005). Under rate-weakening sliding, increases in the
sliding velocity lead to decreases in the basal drag; therefore, a glacier
whose weight is balanced only by basal friction could accelerate unboundedly.
For this reason, rate-weakening sliding has been considered a potential cause
of glacier surges (Lliboutry, 1979; Fowler, 1987). Another implication of
rate-weakening sliding is that the sliding law becomes double-valued, as seen
in Figure 4. This invalidates the commonly used shallow ice approximation of
the Stokes equations which requires the friction law to be invertible (Schoof
& Hewitt, 2013).

We find that, when a steady state along the downsloping section of the curve
is perturbed with a time-dependent water pressure, the cavity quickly evolves
towards the steady state along the upsloping section for a similar value
$\tau_{b,0}/N_{0}$. This is depicted in the lower panel of Figure 5, where we
can see that the volume of the cavity and the sliding speed quickly decrease
and settle around the other steady state associated to $\tau_{b,0}/N_{0}$. We
observed that this phenomenon continues to occur for lower amplitudes and
different frequencies in the oscillations of the effective pressure.
Therefore, our numerical results indicate that under finite perturbations of
the effective pressure, the ice evolves away from the steady states along the
rate-weakening section of the curve. This finding could offer an additional
reason to not use a sliding law with a rate-weakening regime: since such a
regime is unstable in the sense described above, we could expect it to be
unachievable under natural conditions.

Figure 6: Unsteady cavitation around the steady point with
$u_{b,0}/(2ALN_{0}^{n})=0.0761$ with scaled frequencies set to $f=0.4$ (top)
and $f=2$ (lower). On the left panels, we present the evolution of the sliding
speed $u_{b}$ and oscillating effective pressure $N$ (left axis) and of the
cavity volume (right panel). On the right panels, the sliding speed against
the effective pressure; in black, the steady sliding law.

The non-dimensional time it takes a fluid particle to travel across the domain
is of order $t\approx 1$. Therefore, the scaled frequency $f=0.4$ can be
considered a relatively slow frequency that allows the cavity to approximately
follow the steady shapes associated to the effective pressure at each instant
in time as calculated in Section 3. In Figure 6, we compare results obtained
with frequencies $f=0.4$ and $f=2$ to examine the effect of faster
oscillations in the water pressure. With a higher frequency, the magnitude of
the change of cavity volume is significantly lower. Contrastingly, the
amplitude of the sliding speed increases slightly with a higher frequency. The
phase difference between the velocity and the cavity volume also changes when
the frequency is increased. For a high frequency, the maximum velocity is no
longer attained when the cavity reaches its largest extent, but before, when
the cavity is still growing.

Following Sugiyama & Gudmundsson (2004), in Figure 6 we also plot the sliding
speed against the effective pressure throughout one cycle. For the lower
frequency $f=0.4$, a clear loop arises in which the sliding speed is greater
during cavity growth. For the higher frequency of $f=2$, the loop nearly
collapses into a single line. These plots can be compared with those obtained
from field measurements in Sugiyama & Gudmundsson (2004, Figure 5); a
qualitative similarity between both is that higher speeds are reached when $N$
decreases. Additionally, similar plots are presented in Andrews et al. (2014,
Extended Data Figure 4).

##  5 Discussion and conclusions

In this work, we have presented a novel numerical method for solving viscous
contact problems by formulating the Stokes variational inequality as a mixed
problem with Lagrange multipliers. We have also introduced a numerical scheme
for the advection equation specifically designed to ensure both stability and
coherence with the contact conditions of the Stokes problem. This combination
leads to a robust and simple method for solving viscous contact problems. In
Section 3 we have used this method to compute steady cavity configurations. In
particular, we have validated the numerical method by comparing our results
over bedrocks with small amplitudes to the linearised approach in Fowler
(1986) and Schoof (2005) and we have reconstructed steady sliding laws for
different values of the parameter $n$ in Glen’s law (2). Finally, in Section
4, we have explored the temporal evolution of cavities under unsteady
effective pressures and its effect on glacial sliding.

Although not included in this paper, we have compared our method with those
from Gagliardini et al. (2007) and Stubblefield et al. (2021). When computing
the points along the steady sliding law in Figure 4 with our method, the
number of time steps required to converge to a steady state can become very
large (of order 1000) in the highly cavitated stages along the downsloping
region of the curve. This is due to very small scale oscillations that travel
across the cavity but have a significant effect on the calculated values of
$\tau_{b}$ due to the stress singularity at the reattachment point.
Contrastingly, when using the method from Gagliardini et al. (2007), these
oscillations seem to dampen and the method can converge in about 100
iterations for highly cavitated states. We speculate that this is due to the
use of a numerical stabilisation used in Elmer when solving the advection
equation, see (11) in Gagliardini & Zwinger (2008). Despite this difference in
computational times, we find that the basal stress computations carried out
with our method appear to be more accurate due to the exact enforcement of
discrete contact conditions (see Figure 3 and compare with Figure 1 from
Gagliardini et al. (2007)).

We also attempted to solve the subglacial cavity problem with the method from
Stubblefield et al. (2021). Although this penalisation method appears to be
suitable for solving the grounding line and subglacial lake problems presented
in that reference, we found that this method was not able to evolve the cavity
correctly and total detachment would occur after a few time steps. In this
method, via a penalisation term, the constraint
$\boldsymbol{u}\cdot\boldsymbol{n}\leq 0$ in (11) is enforced pointwise at the
discrete level on the attached region. However, due to the discretisation of
the advection equation used in this method, enforcing this condition does not
guarantee that the term $\boldsymbol{u}\cdot\left(\partial\theta/\partial
x,-1\right)$ is exactly zero where $\boldsymbol{u}\cdot\boldsymbol{n}=0$ at
the discrete level. It is this lack of “compatibility” between the
discretisations of the Stokes contact problem and the advection equation that
causes problems when solving the subglacial cavity problem, where the geometry
of the cavity can undergo rapid changes under strong variations of the
effective pressure.

In the numerical method introduced in this paper, we are essentially solving a
Stokes variational inequality by enforcing the average values of
$\boldsymbol{u}_{h}\cdot\boldsymbol{n}$ to be less than or equal to 0 along
the attached region. When formulated as a mixed problem, we arrive at system
(18). However, many other possibilities exist at the discrete level; for
example, we could enforce that the projection of
$\boldsymbol{u}_{h}\cdot\boldsymbol{n}$ onto a space of piecewise linear
functions along the attached region of the boundary be less than or equal to
zero. Alternatively, we could also enforce this constraint to hold for the
values of $\boldsymbol{u}_{h}\cdot\boldsymbol{n}$ at the midpoints of the
edges. We are currently exploring these different discretisations in
combination with different time stepping schemes for the advection equation
(10). We believe that suitable combinations of these methods could yield
further improvements in terms of the speed and robustness of the algorithm.

[Acknowledgements]We acknowledge helpful discussions with O. Gagliardini.

[Funding] This work was supported by the Engineering and Physical Sciences
Research Council (P. E. F., grant numbers EP/V001493/1 and EP/R029423/1). GGdD
was supported by the University of Oxford Mathematical Institute Graduate
Scholarship.

[Declaration of interests]The authors report no conflict of interest.

[Data availability statement]The code for the viscous contact solver presented
in Algorithm 1 is openly available in the archived repository de Diego (2021)
and in
https://bitbucket.org/gonzalogddiego/subglacialcavitysolver_2021/src/master/.
The solver is implemented in Firedrake Rathgeber et al. (2016), using the
version available at zenodo/Firedrake-20211103.0 (2021).

[Author ORCID]G. G. de Diego https://orcid.org/0000-0002-9896-024X; P. E.
Farrell, https://orcid.org/0000-0002-1241-7060; I. J. Hewitt
https://orcid.org/0000-0002-9167-6481

## References

  * Andrews et al. (2014) Andrews, L. C., Catania, G. A., Hoffman, M. J., Gulley, J. D., Lüthi, M. P., Ryser, C., Hawley, R. L. & Neumann, T. A. 2014 Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet. Nature 514 (7520), 80–83.
  * Boffi et al. (2013) Boffi, D., Brezzi, F. & Fortin, M. 2013 Mixed Finite Element Methods and Applications. Springer.
  * de Diego (2021) de Diego, G. G. 2021 Finite element solver for subglacial cavities, https://doi.org/10.5281/zenodo.5647206.
  * de Diego et al. (2021) de Diego, G. G., Farrell, P. E. & Hewitt, I. J. 2021 On the finite element approximation of a semicoercive Stokes variational inequality arising in glaciology, arXiv: 2108.00046.
  * Fowler (1981) Fowler, A. C. 1981 A theoretical treatment of the sliding of glaciers in the absence of cavitation. Philos. T. R. Soc. A 298 (1445), 637–681.
  * Fowler (1986) Fowler, A. C. 1986 A sliding law for glaciers of constant viscosity in the presence of subglacial cavitation. P. Roy. Soc. A-Math. Phy. 407 (1832), 147–170.
  * Fowler (1987) Fowler, A. C. 1987 Sliding with cavity formation. J. Glaciol. 33 (115), 255–267.
  * Gagliardini et al. (2007) Gagliardini, O., Cohen, D., Raback, P. & Zwinger, T. 2007 Finite-element modeling of subglacial cavities and related friction law. J. Geophys. Res.-Earth 112 (F2).
  * Gagliardini & Zwinger (2008) Gagliardini, O. & Zwinger, T. 2008 The ISMIP-HOM benchmark experiments performed using the finite-element code Elmer. The Cryosphere 2 (1), 67–76.
  * Glen (1958) Glen, J. W. 1958 The flow law of ice: a discussion of the assumptions made in glacier theory, their experimental foundation and consequences. In Physics of the Movement of Ice: Symposium at Chamonix 1958, pp. 171–183. Wallingford, UK: Int. Assoc. Hydrol. Sci.
  * Glowinski et al. (1981) Glowinski, R., Lions, J. L. & Trémoliéres, R. 1981 Numerical Analysis of Variational Inequalities. North Holland.
  * Gudmundsson (1997) Gudmundsson, G. H. 1997 Basal-flow characteristics of a non-linear flow sliding frictionless over strongly undulating bedrock. J. Glaciol. 43 (143), 80–89.
  * Helanow et al. (2020) Helanow, C., Iverson, N. R., Zoet, L. K. & Gagliardini, O. 2020 Sliding relations for glacier slip with cavities over three-dimensional beds. Geophys. Res. Lett. 47 (3).
  * Hintermüller et al. (2002) Hintermüller, M., Ito, K. & Kunisch, K. 2002 The primal-dual active set strategy as a semismooth Newton method. SIAM J. Optimiz. 13 (3), 865–888.
  * Hoffman et al. (2016) Hoffman, M. J., Andrews, L. C., Price, S. F., Catania, G. A., Neumann, T. A., Lüthi, M. P., Gulley, J., Ryser, C., Hawley, R. L. & Morriss, B. 2016 Greenland subglacial drainage evolution regulated by weakly connected regions of the bed. Nat. Commun. 7 (1), 13903.
  * Iken (1981) Iken, A. 1981 The effect of the subglacial water pressure on the sliding velocity of a glacier in an idealized numerical model. J. Glaciol. 27 (97), 407–421.
  * Iken & Bindschadler (1986) Iken, A. & Bindschadler, R. A. 1986 Combined measurements of subglacial water pressure and surface velocity of Findelengletscher, Switzerland: Conclusions about drainage system and sliding mechanism. J. Glaciol. 32 (110), 101–119.
  * Kamb (1970) Kamb, B. 1970 Sliding motion of glaciers: Theory and observation. Rev. Geophys. 8 (4), 673–728.
  * LeVeque (2007) LeVeque, R. J. 2007 Finite difference methods for ordinary and partial differential equations: Steady-state and time-dependent problems. SIAM.
  * Lliboutry (1979) Lliboutry, L. 1979 Local friction laws for glaciers: A critical review and new openings. J. Glaciol. 23 (89), 67–95.
  * Mair et al. (2002) Mair, D. W. F., Sharp, M. J. & Willis, I. C. 2002 Evidence for basal cavity opening from analysis of surface uplift during a high-velocity event: Haut Glacier d’Arolla, Switzerland. J. Glaciol. 48 (161), 208–216.
  * Rathgeber et al. (2016) Rathgeber, F., Ham, D. A., Mitchell, L., Lange, M., Luporini, F., McRae, A. T. T., Bercea, G.-T., Markall, G. R. & Kelly, P. H. J. 2016 Firedrake: automating the finite element method by composing abstractions. ACM Trans. Math. Softw. 43 (3), 1–27.
  * Ritz et al. (2015) Ritz, C., Edwards, T. L., Durand, G., Payne, A. J., Peyaud, V. & Hindmarsh, R. C. A. 2015 Potential sea-level rise from Antarctic ice-sheet instability constrained by observations. Nature 528 (7580), 115–118.
  * Schoof (2005) Schoof, C. 2005 The effect of cavitation on glacier sliding. P. Roy. Soc. A-Math. Phy. 461 (2055), 609–627.
  * Schoof (2007) Schoof, C. 2007 Ice sheet grounding line dynamics: Steady states, stability, and hysteresis. J Geophys. Res.-Earth 112 (F3).
  * Schoof (2011) Schoof, C. 2011 Marine ice sheet dynamics. Part 2. a Stokes flow contact problem. J. Fluid Mech. 679, 122–155.
  * Schoof & Hewitt (2013) Schoof, C. & Hewitt, I. J. 2013 Ice-sheet dynamics. Annu. Rev. Fluid Mech. 45 (1), 217–239.
  * Stubblefield et al. (2021) Stubblefield, A. G., Spiegelman, M. & Creyts, T. T. 2021 Variational formulation of marine ice-sheet and subglacial-lake grounding-line dynamics. J. Fluid Mech. 919, A23.
  * Sugiyama & Gudmundsson (2004) Sugiyama, S. & Gudmundsson, G. H. 2004 Short-term variations in glacier flow controlled by subglacial water pressure at Lauteraargletscher, Bernese Alps, Switzerland. J. Glaciol. 50 (170), 353–362.
  * Weertman (1957) Weertman, J. 1957 On the sliding of glaciers. J. Glaciol. 3 (21), 33–38.
  * zenodo/Firedrake-20211103.0 (2021) zenodo/Firedrake-20211103.0 2021 Software used in ‘Numerical approximation of viscous contact problems applied to glacial sliding’, https://doi.org/10.5281/zenodo.5643646.
  * Zoet & Iverson (2015) Zoet, L. K. & Iverson, N. R. 2015 Experimental determination of a double-valued drag relationship for glacier sliding. J. Glaciol. 61 (225), 1–7.