# COMET-QE and Active Learning for Low-Resource Machine Translation

Everlyn Asiko Chimoto
University of Cape Town, South Africa
African Institute for Mathematical
Sciences
<EMAIL_ADDRESS>
Bruce A. Bassett
University of Cape Town, South Africa
African Institute for Mathematical
Sciences, South Africa
South African Astronomical Observatory
<EMAIL_ADDRESS>

###### Abstract

Active learning aims to deliver maximum benefit when resources are scarce. We
use COMET-QE, a reference-free evaluation metric, to select sentences for low-
resource neural machine translation. Using Swahili, Kinyarwanda and Spanish
for our experiments, we show that COMET-QE significantly outperforms two
variants of Round Trip Translation Likelihood (RTTL) and random sentence
selection by up to 5 BLEU points for 20k sentences selected by Active Learning
on a 30k baseline. This suggests that COMET-QE is a powerful tool for sentence
selection in the very low-resource limit.

##  1 Introduction

Active Learning (AL) is a technique that involves querying an Oracle for
labels based on a selection strategy. For Neural Machine Translation (NMT),
active learning involves selecting sentences for translation based on
estimated model uncertainty on the sentence or the diversity of the sentence.
The aim of AL is to choose sentences that would maximally improve the NMT
model. This is particularly attractive for low-resource languages where
budgets are small, which includes most of Africa’s languages which have
historically been neglected Nekoto et al. (2020).

There have been several advances in Active Learning for NMT with Round Trip
Translation Likelihood (RTTL) being the current state of the art (Zeng et al.,
2019; Peris and Casacuberta, 2018; Zhao et al., 2020; Koneru et al., 2022).
These advancements have typically been tested by artificially restricting
high-resource languages which mimic true low-resource languages. There have
been few empirical evaluations of how these AL techniques perform on a low-
resource African language with less than 100K sentences. We tackle this
problem while simultaneously proposing a new AL sentence selection strategy
based on COMET-QE (Rei et al., 2020) and show that it outperforms two versions
of RTTL, one of the most successful AL selection algorithms, on Swahili,
Spanish and Kinyarwanda. Interestingly the RTTL variants perform worse than
random sentence sampling in this context.

##  2 Active Learning for NMT

There are several AL sentence selection frameworks including pool-based AL,
membership query synthesis and stream-based AL Settles (2009). We focus on
pool-based AL where we assume their exist a large pool of monolingual
sentences for the language. Pool-based AL in NMT is an iterative process that
involves selecting batches of sentences after each given AL cycle (Haffari et
al., 2009; Bloodgood and Callison-Burch, 2010). Given a set of parallel
sentences $\mathcal{L}$, a model $\mathcal{M}$ is trained. $\mathcal{M}$ is
then used to translate monolingual data, $\mathcal{U}$ and a selection
strategy scores and ranks the sentences in $\mathcal{U}$. The top sentences
are then selected for translation. After translation, the selected sentences
are added to set $\mathcal{L}$ and used to retrain and improve the model
(Ambati et al., 2010); see Figure 1. This process is repeated until a
specified stopping criteria is fulfilled, e.g. total translation budget
reached.

Figure 1: Active Learning framework where the labelled set $\mathcal{L}$ is
used to train model $\mathcal{M}$ which score monolingual sentences based on
an uncertainty measure. After scoring, the sentences selected based on the
annotation budget get translated and added to the training set for retraining.

The selection strategy aims to choose the N sentences that will improve the
model the most. Selection strategies are typically either model-driven or
data-driven. In Model-driven techniques a model scores sentences whereas Data-
driven methods, also known as model-agnostic, involve comparing the set of
labelled and unlabelled vernacular sentences to score the sentences. The
current State-of-the-art for model-driven selection is RTTL Zeng et al.
(2019); Koneru et al. (2022) which we discuss now.

###  2.1 Round Trip Translation Likelihood (RTTL)

Zeng et al. (2019) proposed RTTL as an AL selection strategy. RTTL involves
training 2 models: i) Model $\mathcal{M}$ that translates from the source to
target language and ii) Model $\mathcal{M}_{rev}$ that reverse translates back
to the source language. Given a source sentence $x$, $\mathcal{M}$ generates a
translation $\hat{y}$. $\mathcal{M}_{rev}$ translates $\hat{y}$ back to the
source language. However, instead of comparing of $x$ and $\hat{x}$, the
probability of getting the actual sentence $x$ can be used to select sentences
for translation. In RTTL, the key likelihood is defined in terms of $x$ and
$\hat{y}$ as:

$\phi(x,\mathcal{M},\mathcal{M}_{rev})=-\frac{1}{L}\log P_{rev}(x|\hat{y})$

where $L$ is the length of the sentence and $P_{rev}$ is the probability
generated by the reverse model which estimates the probability of the original
sentence given the translated sentence. The lower the probability the more
uncertain the model is in translating the sentence and the more likely it will
thus be selected for translation. We refer to this version as RTTL and our
variant, discussed below.

###  2.2 Random Baseline

In addition to RTTL we compare our new sentence selection methods against a
random baseline in which sentences are selected at random.

##  3 New Active Learning Algorithms

In this section, we detail our proposed new sentence selection techniques.

###  3.1 COMET-QE

Estimating the quality of a translation is a key area in evaluating machine
translation. There have been several advancements in estimating the quality of
a translation without having the reference or label. These reference-free
metrics can be used as a selection strategy. Just as with RTTL, if the quality
of a model translation is poor, this sentence is a good candidate for
selection by active learning.

COMET-QE is a reference-free metric proposed by Rei et al. (2020) under the
WMT 2020 Shared Task on Metrics. COMET-QE correlates well to human assessments
of translation quality.

COMET-QE uses a pre-trained model, XLM-RoBERTa (XLM-R) (Conneau et al., 2020)
to encode both the source sentences as well as the translation. Thereafter, a
pooling layer converts the embeddings into segment-level vectors and finally
the last layer produces a quality score; see Figure 2.

For our current work, Swahili, Spanish and Kinyarwanda source sentences and
English hypotheses were fed into the COMET-QE model, generating quality
estimates for each pair of sentences. This AL algorithm then selects the
lowest quality candidate sentence pairs for translation by the annotator.

Figure 2: The COMET-QE architecture. The model receives a source sentence, s
and a hypothesis, h, passes them through the pre-trained XLM-Roberta model and
performs regression to estimate a quality score. Figure source: Unbabel .

###  3.2 Stratified RTTL

Since RTTL selects sentences with low likelihood given the reverse translation
model, one may be concerned that it selects overly long/short sentences that
do not truly add diversity to the training dataset. Thus we propose stratified
RTTL (S-RTTL) which implements RTTL on length-based stratification of the test
or training sets.

S-RTTL looks at the distribution of sentence length in the test or training
set and selects the highest scoring sentences in each length bin/stratum
proportional to the number of test/training sentences with that length. This
ensures that RTTL is not biased to the sentences of a certain length. In
effect, given the initial training data distribution, S-RTTL would select
sentences that fall under the same length distribution of test/training. An
example, if the training data has 10% sentences of length 1-10, S-RTTL picks
the most uncertain sentences of length 1-10 to fill 10% of sentences it
selects.

When matching the test distribution (if it is available) S-RTTL provides a way
to help deal with dataset shift. Here we apply S-RTTL to match the training
data length distribution, though because our test/training split is performed
randomly there is no functional difference.

##  4 Experiments

###  4.1 Datasets

We run experiments translating three languages into English, namely
Kinyarwanda, Spanish and Swahili. We choose Swahili and Kinyarwanda because
they are true low-resource African languages whereas we chose Spanish because
it was the original language used to demonstrate RTTL, though we mimic a very
low-resource setting in our case.

Our experiments on Swahili-English uses the SAWA Swahili dataset which
consists of data from various domains (De Pauw et al., 2009; Pauw et al.,
2011). The SAWA dataset is provided on request by the authors (De Pauw et al.,
2009). In total, there were 272,544 parallel sentences. Pre-processing
included: removing missing rows with missing translations and removing
sentences with more than 100 words. This clean set consisted of 89,505
parallel sentences. We carried out 5 runs where we split the data into 5
separate folds with each test consisting of 17,901 parallel sentences. Of the
remaining sentences, we select 30,000 sentences as our base training set and
1,000 sentences as our validations set. This leaves 40,604 sentences to select
from using AL. With each iteration, we use AL to select 5,000 new sentences
and add them to the training set. The validation and test set remain the same
for each active learning strategy through all iterations to ensure fair
comparison of the methods. Byte Pair Encoding with 4,000 merge operations was
applied to the dataset (Sennrich et al., 2016).

For the Spanish-English dataset, we utilised Europarl version 7 which consists
of $\sim 2M$ parallel sentences. For the Kinyarwanda-English dataset, we used
the JW300 dataset (Agić and Vulić, 2019) which contains 436,753 parallel
sentences. Pre-processing for these 2 datasets included removing sentences
with more than 100 words and removing sentences that are the same in both
languages. We then randomly chose 148,000 sentences to work with. As with the
Swahili-English dataset, we split the data into 5 separate folds. With each
test fold, we randomly selected 17,000 sentences as the test set. Thereafter,
we concatenate the remaining sentences with the other sentences in the
training fold. We randomly select 30,000 sentences as our base training set
and 1,000 sentences as our validations set. Contrary to the Swahili-English
dataset, due to computational limitations we only conduct one AL iteration
where we select 20k sentence out of the 100k remaining pool of training
sentences.

###  4.2 Training

We used the transformer-based NMT model JoeyNMT library (Kreutzer et al.,
2019). The configuration for all the models were: 6 layers for both encoder
and decoder, 4 heads, embedding dimension of 256, feedforward size of 1,024.
We train over 60 epochs, applying Adam optimizer with $\beta_{1}$ =0.9 and
$\beta_{2}$=0.98 and a learning rate of 0.001.

All training was done on Google Colab which typically provides the Tesla K80
GPU. Every model iteration, both training and testing, was run an average of 4
hours with approximately 12M parameters.

Figure 3: Comparison of the four AL approaches used on Swahili in terms of
$\Delta$ BLEU score relative to the baseline of 30,000 sentences. COMET-QE as
an AL strategy outperforms the rest of the techniques whereas the RTTL
performs the poorest.

###  4.3 Results and Discussion

| Sentence Length | Number of symbols
---|---|---
Iteration | Random | RTTL | S- RTTL | COMET-QE | Random | RTTL | S-RTTL | COMET-QE
35k | 11 | 5.3 | 10.2 | 24.7 | 12.4 | 5.5 | 11.3 | 29.0
40k | 10.7 | 9.0 | 10.4 | 16.4 | 12.1 | 9.9 | 11.5 | 18.7
45k | 10.7 | 9.8 | 10.2 | 12.4 | 12.0 | 10.8 | 11.2 | 14.0
50k | 10.8 | 11.0 | 10.5 | 9.8 | 12.1 | 12.3 | 11.7 | 10.9
55k | 10.9 | 12.7 | 10.8 | 7.7 | 12.3 | 14.5 | 12.1 | 8.4
Average | 10.82 | 9.56 | 10.42 | 14.2 | 12.18 | 10.6 | 11.56 | 16.2

Table 1: Analysis of the length of sentences and number of symbols in a
sentence picked at each iteration of the AL process for Swahili. COMET-QE
favors longer sentences while RTTL favors shorter sentences. RTTL seems to
also pick sentences with a high ratio between sentence length and number of
symbols.

Despite its known limitations we choose the BLEU score (Papineni et al., 2002)
as a standard metric for comparing the results for the different models. We
find that COMET-QE outperforms all of the other methods (RTTL, S-RTTL and
random) on all three languages, often by a significant margin. This is shown
for Swahili in Figure (3) and for all languages in Table (2). COMET-QE has
nearly double the BLEU score gain compared to RTTL for Swahili after adding
25,000 sentences to the baseline. Error contours are estimated 1-$\sigma$
bounds based on our multiple runs.

Our results show that COMET-QE is a promising active learning selection
strategy in the very low-resource setting. Looking deeper into the sentences
that AL with COMET-QE selects (see Table 1), we see that COMET-QE initially
selects sentences that are significantly longer (24.7 words per sentence) with
more symbols than any of the other methods. With each iteration it
progressively chooses shorter and shorter sentences, ending with just 7.7
words per sentence.

The exact opposite strategy is chosen by RTTL which starts off selecting short
sentences (5.3 words/sentence) and gradually lengthening them, ending with
12.7 words per sentence. Random and S-RTTL choose a middle path, with the
average length and symbol complexity being roughly constant with each
iteration (about 10.5 words/sentence).

This same pattern is seen in the number of unique words in each group of 5,000
sentences chosen by the algorithms at each iteration. COMET-QE’s initial
choice has almost 20,000 unique words, nearly double that of RTTL. In
contrast, by the 5th iteration, COMET-QE has halved the number of unique words
while RTTL has increased by $\sim 20\%$, see Table (3).

Language | Random | RTTL | S-RTTL | COMET-QE
---|---|---|---|---
Swahili | $5.7\pm 0.8$ | $3.0\pm 1.2$ | $4.1\pm 0.9$ | $\bf{7.7\pm 1.2}$
Spanish | $3.2\pm 0.4$ | $2.8\pm 0.3$ | $1.7\pm 0.5$ | $\bf{3.5\pm 0.4}$
Kinyarwanda | $3.4\pm 0.3$ | $2.3\pm 0.3$ | $1.5\pm 0.5$ | $\bf{5.1\pm 0.5}$

Table 2: Gains in BLEU score relative to a 30k sentence baseline for each of
Swahili, Spanish and Kinyarwanda with 20k sentences selected using the various
AL algorithms. Error estimates are $1-\sigma$ errors over the five runs. The
best algorithm is highlighted in bold: in all cases COMET-QE outperforms the
other AL techniques. With these small amounts of data RTTL and S-RTTL perform
the poorest. The baseline scores for the languages with 30k sentences were:
Swahili: $20.8\pm 0.7$, Spanish: $29.0\pm 0.2$, Kinyarwanda: $15.7\pm 1.2$.

In Figure (3) we notice that RTTL performs the worst of all the AL strategies,
contrary to what was reported in Zeng et al. (2019) and Koneru et al. (2022)
where RTTL significantly outperforms random. However, poor performance of RTTL
relative to random selection was also recently observed in Hu and Neubig
(2021). Here we attribute the poor performance of RTTL to the fact that we are
exploring the very-low resource setting (only 50k sentences in total). RTTL
requires using the learned translation models in both directions, e.g. En-Sw
and Sw-En, and hence errors in each direction compound and are amplified. In
comparison COMET-QE only requires translation in one direction and therefore
one might expect it to perform better in the very low-resource setting.
Exploring the selections made by RTTL, we note that it tended to select
sentences that are not appropriate because they are dominated by symbols or
other noise of no value to learning a translation model.

Interestingly, S-RTTL improves on RTTL but despite choosing sentences with
almost identical average sentence length to the random selection algorithm,
still returns a BLEU score about 1.5 lower on average compared to random; see
Table (1).

| Random | RTTL | S-RTTL | COMET-QE
---|---|---|---|---
35k | 12,448 | 10,662 | 14,356 | 19,995
40k | 12,253 | 12,797 | 13,375 | 15,755
45k | 12,299 | 12,709 | 12,719 | 13,258
50k | 12,227 | 12,779 | 12,397 | 11,589
55k | 12,394 | 13,095 | 12,107 | 9,948

Table 3: Number of unique Swahili words picked at each iteration by each of
the algorithms.

## Conclusions

We study various active learning strategies in a very low-resource Neural
Machine Translation (NMT) setting, specifically using Swahili, Kinyarwanda and
Spanish to English as our target languages. Our experiments show that using
COMET-QE, a reference-free quality estimation metric, as an Active Learning
strategy significantly performs Round Trip Translation Likelihood (RTTL), the
current SOTA active learning method as well as random selection.

COMET-QE initially selects long sentences with a much greater diversity of
words than the other algorithms, leading to a BLEU score increase of up to 5
over the 30k baseline for Swahili with just an extra 5,000 sentences added.
COMET-QE then progressively chooses shorter sentences with fewer rare words
with each successive iteration. Conversely, we find that RTTL, performs worse
than random sentence selection. RTTL starts by selecting short sentences and
then gradually increases the sentence length in each iteration. To improve its
performance, we introduced Stratified RTTL (S-RTTL) which matches the test or
training set sentence length distribution. This method leads to improved
performance for Swahili but not for Spanish and Kinyarwanda and was still
worse than random. We attribute the poor performance of RTTL and S-RTTL to the
very low-resource setting: the translation models are still weak due to lack
of data and since RTTL uses them in both translation directions its estimates
of sentence quality are weak.

This shows that advances of translation quality estimation metrics can be
effectively leveraged to improve active learning strategies for very low-
resource NMT.

## Limitations

The above experiments were conducted on only three languages and limit the
study to a very low-resource setting, below 55k sentences. This means that the
results may not apply to other languages and investigating the performance of
COMET-QE as well as RTTL on other languages and settings is left to future
work.

## Acknowledgements

We wish to thank the Masakhane community for support. We thank Emmanuel
Dufourq and the anonymous reviewers for comments and contributions to the
draft. This publication was made possible by a grant from Carnegie Corporation
of New York (provided through the African Institute for Mathematical
Sciences). The statements made and views expressed are solely the
responsibility of the authors.

## References

  * Agić and Vulić (2019) Željko Agić and Ivan Vulić. 2019.  JW300: A wide-coverage parallel corpus for low-resource languages.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 3204–3210, Florence, Italy. Association for Computational Linguistics.
  * Ambati et al. (2010) Vamshi Ambati, Stephan Vogel, and Jaime Carbonell. 2010.  Active learning and crowd-sourcing for machine translation.  In _Proceedings of the Seventh International Conference on Language Resources and Evaluation (LREC’10)_ , Valletta, Malta. European Language Resources Association (ELRA).
  * Bloodgood and Callison-Burch (2010) Michael Bloodgood and Chris Callison-Burch. 2010.  Bucking the trend: Large-scale cost-focused active learning for statistical machine translation.  In _Proceedings of the 48th Annual Meeting of the Association for Computational Linguistics_ , pages 854–864, Uppsala, Sweden. Association for Computational Linguistics.
  * Conneau et al. (2020) Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. 2020.  Unsupervised cross-lingual representation learning at scale.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 8440–8451, Online. Association for Computational Linguistics.
  * De Pauw et al. (2009) Guy De Pauw, Peter Waiganjo Wagacha, and Gilles-Maurice de Schryver. 2009.  The SAWA corpus: A parallel corpus English - Swahili.  In _Proceedings of the First Workshop on Language Technologies for African Languages_ , pages 9–16, Athens, Greece. Association for Computational Linguistics.
  * Haffari et al. (2009) Gholamreza Haffari, Maxim Roy, and Anoop Sarkar. 2009.  Active learning for statistical phrase-based machine translation.  In _Proceedings of Human Language Technologies: The 2009 Annual Conference of the North American Chapter of the Association for Computational Linguistics_ , pages 415–423, Boulder, Colorado. Association for Computational Linguistics.
  * Hu and Neubig (2021) Junjie Hu and Graham Neubig. 2021.  Phrase-level active learning for neural machine translation.  In _Proceedings of the Sixth Conference on Machine Translation_ , pages 1087–1099, Online. Association for Computational Linguistics.
  * Koneru et al. (2022) Sai Koneru, Danni Liu, and Jan Niehues. 2022.  Cost-effective training in low-resource neural machine translation.  _CoRR_ , abs/2201.05700.
  * Kreutzer et al. (2019) Julia Kreutzer, Jasmijn Bastings, and Stefan Riezler. 2019.  Joey NMT: A minimalist NMT toolkit for novices.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP): System Demonstrations_ , pages 109–114, Hong Kong, China. Association for Computational Linguistics.
  * Nekoto et al. (2020) Wilhelmina Nekoto, Vukosi Marivate, Tshinondiwa Matsila, Timi Fasubaa, Taiwo Fagbohungbe, Solomon Oluwole Akinola, Shamsuddeen Muhammad, Salomon Kabongo Kabenamualu, Salomey Osei, Freshia Sackey, Rubungo Andre Niyongabo, Ricky Macharm, Perez Ogayo, Orevaoghene Ahia, Musie Meressa Berhe, Mofetoluwa Adeyemi, Masabata Mokgesi-Selinga, Lawrence Okegbemi, Laura Martinus, Kolawole Tajudeen, Kevin Degila, Kelechi Ogueji, Kathleen Siminyu, Julia Kreutzer, Jason Webster, Jamiil Toure Ali, Jade Abbott, Iroro Orife, Ignatius Ezeani, Idris Abdulkadir Dangana, Herman Kamper, Hady Elsahar, Goodness Duru, Ghollah Kioko, Murhabazi Espoir, Elan van Biljon, Daniel Whitenack, Christopher Onyefuluchi, Chris Chinenye Emezue, Bonaventure F. P. Dossou, Blessing Sibanda, Blessing Bassey, Ayodele Olabiyi, Arshath Ramkilowan, Alp Öktem, Adewale Akinfaderin, and Abdallah Bashir. 2020.  Participatory research for low-resourced machine translation: A case study in African languages.  In _Findings of the Association for Computational Linguistics: EMNLP 2020_ , pages 2144–2160, Online. Association for Computational Linguistics.
  * Papineni et al. (2002) Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zhu. 2002.  Bleu: a method for automatic evaluation of machine translation.  In _Proceedings of the 40th Annual Meeting of the Association for Computational Linguistics_ , pages 311–318, Philadelphia, Pennsylvania, USA. Association for Computational Linguistics.
  * Pauw et al. (2011) Guy De Pauw, Peter Waiganjo Wagacha, and Gilles-Maurice de Schryver. 2011.  Exploring the sawa corpus: collection and deployment of a parallel corpus english—swahili.  _Language Resources and Evaluation_ , 45:331–344.
  * Peris and Casacuberta (2018) Álvaro Peris and Francisco Casacuberta. 2018.  Active learning for interactive neural machine translation of data streams.  In _Proceedings of the 22nd Conference on Computational Natural Language Learning_ , pages 151–160, Brussels, Belgium. Association for Computational Linguistics.
  * Rei et al. (2020) Ricardo Rei, Craig Stewart, Ana C Farinha, and Alon Lavie. 2020.  Unbabel’s participation in the WMT20 metrics shared task.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 911–920, Online. Association for Computational Linguistics.
  * Sennrich et al. (2016) Rico Sennrich, Barry Haddow, and Alexandra Birch. 2016.  Neural machine translation of rare words with subword units.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 1715–1725, Berlin, Germany. Association for Computational Linguistics.
  * Settles (2009) Burr Settles. 2009.  Active learning literature survey.
  * (17) Unbabel. Accessed June 2022.  Comet metrics.  Github https://unbabel.github.io/COMET/html/models.html.
  * Zeng et al. (2019) Xiangkai Zeng, Sarthak Garg, Rajen Chatterjee, Udhyakumar Nallasamy, and Matthias Paulik. 2019.  Empirical evaluation of active learning techniques for neural MT.  In _Proceedings of the 2nd Workshop on Deep Learning Approaches for Low-Resource NLP (DeepLo 2019)_ , pages 84–93, Hong Kong, China. Association for Computational Linguistics.
  * Zhao et al. (2020) Yuekai Zhao, Haoran Zhang, Shuchang Zhou, and Zhihua Zhang. 2020.  Active learning approaches to enhancing neural machine translation.  In _Findings of the Association for Computational Linguistics: EMNLP 2020_ , pages 1796–1806, Online. Association for Computational Linguistics.


The effect of model freshwater flux biases

on the multi-stable regime of the AMOC

Henk A. Dijkstra 111Correspondence: Email<EMAIL_ADDRESS>and René M. van
Westen

Institute for Marine and Atmospheric research Utrecht

Department of Physics, Utrecht University

Princetonplein 5, 3584 Utrecht, the Netherlands

Submitted to Tellus A

Version of

###### Abstract

It is known that global climate models (GCMs) have substantial biases in the
surface freshwater flux which forces the ocean component of these models.
Using numerical bifurcation analyses on a global ocean model, we study here
the effect of a specific freshwater flux bias on the multiple equilibrium
regime of the Atlantic Meridional Overturning Circulation (AMOC). We find that
a (positive) freshwater flux bias over the Indian Ocean shifts the multiple
equilibrium regime to larger values of North Atlantic freshwater input but
hardly affects the associated hysteresis width. The magnitude of this shift
depends on the way the anomalous North Atlantic freshwater flux is
compensated. We explain the changes in bifurcation diagrams using the
freshwater balance over the Atlantic basin. The results suggest that state-of-
the-art GCMs may have an AMOC multiple equilibrium regime, but that it is
located in a parameter regime that is considered unrealistic and hence is not
explored.

##  1 Introduction

The AMOC has been proposed as one of the tipping elements in the climate
system Lenton et al. (2008); Armstrong McKay et al. (2022), indicating that it
may undergo a relatively rapid change under a slowly developing forcing. The
AMOC is thought to be particularly sensitive to its freshwater forcing, either
through the surface freshwater flux, through input of freshwater due to ice
melt (e.g. from the Greenland Ice Sheet) or through river runoff Rahmstorf et
al. (2005). These freshwater fluxes affect the ocean density differences which
control the AMOC strength and change the heat and salt transport carried by
AMOC Marotzke (2000).

From conceptual models, the tipping behaviour is clearly related to the multi-
stability properties of the AMOC. For example, in the Stommel two-box model
Stommel (1961), there is an interval of surface freshwater forcing where two
stable steady AMOC states exist and tipping occurs due to transitions between
these states. The central feedback responsible for the tipping behaviour is
the salt-advection feedback, where a freshwater perturbation in the North
Atlantic causes a weakening of the AMOC which leads to less northward salt
transport and hence amplification of the perturbation Marotzke (2000).

Precise boundaries of the multiple-equilibrium regime of the AMOC have been
obtained using conceptual models Cessi (1994); Cimatoribus et al. (2012) and
fully-implicit ocean-climate models De Niet et al. (2007); Toom et al. (2012);
Mulder et al. (2021). One of the important results of these studies is that
the existence of the multiple-equilibrium regime can be related (in these
models) to an observable quantity Rahmstorf (1996), now often called the AMOC
stability (or regime) indicator. This indicator has many different notations
in the literature, e.g. $M_{ov}$ De Vries and Weber (2005) or $F_{ov}$ Hawkins
et al. (2011). Here, we will follow Weijer et al. (2019) and use $F_{ovS}$
($F_{ovN}$) as the freshwater transport carried by the AMOC over the southern
(northern) boundary at 35∘S (60∘N) of the Atlantic basin Dijkstra (2007);
Huisman et al. (2010); Liu et al. (2017). Available observations Bryden et al.
(2011) show that the present-day AMOC is exporting freshwater out of the
Atlantic ($F_{ovS}<0$). It is known that $F_{ovS}$ ignores some relevant
processes Gent (2018), but if one accepts that $F_{ovS}$ is a proper
indicator, the AMOC is in a multiple-equilibrium regime based its observed
values Weijer et al. (2019).

Less precise estimates of the multiple-equilibrium regime boundaries can be
obtained from global climate models. In so-called quasi-equilibrium
experiments, the freshwater forcing is changed very slowly such that the model
state stays close to the (slowly changing) equilibrium. When the freshwater
forcing is varied in both directions and covers the multiple-equilibrium
regime, regime boundaries can be inferred from the so-called hysteresis width,
i.e. the freshwater forcing values where the AMOC collapses and recovers. The
rate of forcing is important here and if this is much faster than the
equilibration time scale of the steady state, the approximations of the regime
boundaries become worse and also rate-induced tipping may occur Lohmann et al.
(2021). Such quasi-equilibrium experiments have been performed with many Earth
System Models of Intermediate Complexity (EMICs) Rahmstorf et al. (2005) and
the FAMOUS model, the latter being a GCM with a relatively coarse horizontal
resolution of $2.5^{\circ}\times 3.75^{\circ}$ Hawkins et al. (2011). The
Community Climate System Model (CCSM3) is probably the most detailed model in
which AMOC hysteresis behavior has been investigated Hu et al. (2012).

Mostly due to computational constraints, often only the AMOC response to
particular freshwater forcing perturbations is considered in state-of-the-art
GCMs. In these so-called ‘hosing experiments’ Stouffer et al. (2006), quite a
diversity of model behavior is found. It is not known whether a multiple-
equilibrium AMOC regime exists in such models and only sporadic indications of
such a regime have been found Mecking et al. (2016); Jackson and Wood (2018a,
b). The problem is that it is difficult to assess whether the weak AMOC states
computed are equilibrium solutions of the models. What these model studies
certainly have shown is that an AMOC weakening would have severe impacts on
the climate system, affecting sea level and regional temperatures in many
areas around the world Vellinga et al. (2002); Jackson et al. (2015); Liu et
al. (2020); Orihuela-Pinto et al. (2022).

Of course, the real present-day AMOC may not have a multiple equilibrium
regime and GCMs may model that correctly. On the other hand, the real AMOC may
be in such a regime, and the GCMs may not capture it. In that case, the GCMs
misrepresent (or miss) crucial processes, such that they do not display
tipping behaviour. A prominent example is the incorrect representation of
ocean eddy transport processes in GCMs, which may prevent the existence of a
multi-stable AMOC regime. Another possibility is that GCMs capture the
relevant processes but the parameters in the models are not correct, such that
a multi-stable regime does not occur due to model biases. An example of this
is that many GCMs are considered to have a too stable AMOC due to biases in
the freshwater transport in the Atlantic Ocean Drijfhout et al. (2013);
Mecking et al. (2016).

In van Westen and Dijkstra (2023), it was shown that several persistent biases
in state-of-the art climate models, of the Coupled Model Intercomparison
Projects (CMIP) version 6, lead to an AMOC with an Atlantic freshwater
transport that is in disagreement with observations (i.e., many models have
$F_{ovS}>0$). The most important model bias is the surface freshwater flux
over the Indian Ocean, which affects the freshwater transport at 34∘S in the
Atlantic through Agulhas Leakage. Motivated by these results, we revisit the
bifurcation analysis on the global ocean-climate model Dijkstra (2007) and
determine the effects of Indian Ocean surface freshwater flux biases on the
multi-stable regime of the AMOC. In section 2, the fully-implicit model used
is shortly summarised and the continuation methodology to compute bifurcation
diagrams is described. Then in section 3, we focus on the effect of freshwater
biases on the bifurcation diagrams. Mechanisms of the shift in the multi-
stable regimes are analysed using the freshwater balance over the Atlantic. A
summary and discussion follows in section 4.

##  2 Formulation

###  2.1 Model

The fully-implicit global ocean model used in this study is described in
detail in Dijkstra and Weijer (2005) to which the reader is referred for full
details. The governing equations of the ocean model are the hydrostatic,
primitive equations in spherical coordinates on a global domain which includes
continental geometry as well as bottom topography. The ocean velocities in
eastward (zonal) and northward (meridional) directions are indicated by $u$
and $v$, the vertical velocity is indicated by $w$, the pressure by $p$ and
the temperature and salinity by $T$ and $S$, respectively. The horizontal
resolution of the model is about $4^{\circ}$ (a $96\times 38$ Arakawa C-grid
on a domain [$180^{\circ}$ W$,180^{\circ}$E] $\times$ [$85.5^{\circ}$
S$,85.5^{\circ}$N]) and the grid has 12 vertical levels. The vertical grid is
non-equidistant with the surface (bottom) layer having a thickness of 50 m
(1000 m), respectively.

Vertical and horizontal mixing of momentum and of tracers (i.e., heat and
salt) are represented by a Laplacian formulation with prescribed ‘eddy’
viscosities $A_{H}$ and $A_{V}$ and eddy diffusivities $K_{H}$ and $K_{V}$,
respectively. As in Dijkstra (2007), we will use the depth dependent values of
$K_{V}$ and $K_{H}$ Bryan and Lewis (1979); England (1993) given by

$K_{V}(z)&=K^{0}_{V}-A_{s}\arctan(\lambda_{V}(z-z_{*})),\\\
K_{H}(z)=K^{0}_{H}+(A_{r}-K^{0}_{H})e^{\frac{z}{\lambda_{H}}},$ (1)

with $z\in[-5000,0]$ m. Here, $K^{0}_{H}=0.5\times 10^{3}$ m2s-1,
$A_{r}=1.0\times 10^{3}$ m2s-1, $K^{0}_{V}=8.0\times 10^{-5}$ m2s-1,
$A_{s}=3.3\times 10^{-5}$ m2s-1, $\lambda_{V}=4.5\times 10^{-3}$ m-1,
$\lambda_{H}=5\times 10^{2}$ m and $z_{*}=-2.5\times 10^{3}$ m. A plot of the
vertical structure of $K_{V}$ and $K_{H}$ can be found in Figure 1 of Dijkstra
(2007). In this way, the vertical diffusivity $K_{V}$ increases from
$0.31\times 10^{-4}$ m2s-1 at the surface to $1.3\times 10^{-4}$ m2s-1 near
the bottom of the flow domain. The horizontal diffusivity $K_{H}$ increases
monotonically from $0.5\times 10^{3}$ m2s-1 at the bottom of the ocean to
$1.0\times 10^{3}$ m2s-1 near the surface.

The ocean flow is forced by the observed annual-mean wind stress as given in
Trenberth et al. (1989). The upper ocean is coupled to a simple energy-balance
atmospheric model (see Appendix in Dijkstra and Weijer (2005)) in which only
the heat transport is modelled (no moisture transport). The freshwater flux
will be prescribed in each of the results in section 3 and the model has no
sea-ice component. The surface forcing is represented as a body forcing over
the upper layer. On the continental boundaries, no-slip conditions are
prescribed and the heat- and salt fluxes are zero. At the bottom of the ocean,
both the heat and salt fluxes vanish and slip conditions are assumed.

###  2.2 Methods

The discretised steady equations can be written as a nonlinear algebraic
system of equations of the form

${\bf G}({\bf x},\mu)=0,$ (2)

where ${\bf x}$ is the state vector and $\mu$ is one of the parameters of the
model. For the global ocean model (with a 4∘ horizontal resolution and 12
layers in the vertical) the dimension of the state space (and of $\bf x$) is
$96\times 38\times 13\times 6=284,544$; where the number $13$ comes from the
$12$ ocean levels plus the atmospheric energy balance model.

We use pseudo-arclength continuation Keller (1977), where the branch of steady
solutions versus $\mu$ is parametrised by an arclength $s$. To close the set
of equations (because of the new variable $s$) the arclength is normalised
leading to the equations

${\bf G}({\bf x}(s),\mu(s))&=0,\\\ \dot{\bf x}^{T}_{0}({\bf x}(s)-{\bf
x}_{0})+\dot{\mu}_{0}(\mu(s)-\mu_{0})-(s-s_{0})=0,$ (3)

where $({\bf x}_{0},\mu_{0})$ is a previously computed solution and the dot
indicates differentiation to $s$. The linear stability of each steady state is
determined by solving a generalized eigenvalue problem using the Jacobi-
Davidson QZ method Dijkstra (2005).

Figure 1: _Areas where freshwater flux anomalies are applied with their
strengths $\gamma_{A}$ and $\gamma_{I}$; also the global compensation region
(a) and the Pacific region (b) is shown. _

The procedure to compute bifurcation diagrams of the model, including biases
in the freshwater forcing, is the following:

  * (i)

Under restoring conditions for the surface salinity field Levitus (1994), a
steady solution is determined for standard values of the parameters of the
model Dijkstra and Weijer (2005). From this steady solution the freshwater
flux, below referred to the Levitus flux $F^{L}_{S}$, is diagnosed.

  * (ii)

A freshwater flux over a region near New Foundland (Fig. 1) with domain
[$60^{\circ}$ W$,24^{\circ}$ W] $\times$ [$54^{\circ}$ N $,66^{\circ}$ N] is
prescribed (in addition to $F^{L}_{S}$) with strength $\gamma_{A}F^{A}_{S}$
Sv, where $F^{A}_{S}=1$ in this domain and zero outside. Similarly, a bias
freshwater flux is prescribed over the Indian Ocean domain [$52^{\circ}$
E$,104^{\circ}$ E] $\times$ [$20^{\circ}$ S, $10^{\circ}$ N] with amplitude
$\gamma_{I}F^{I}_{S}$ Sv, where $F^{I}_{S}=1$ in this domain and zero outside.
The total freshwater flux is prescribed as

$F_{S}=F^{L}_{S}+\gamma_{A}F^{A}_{S}+\gamma_{I}F^{I}_{S}-QF^{C}_{S},$ (4)

where $F^{C}_{S}=1$ in a compensation domain (specified below) and the
quantity $Q$ is determined such that

$\int_{S_{oa}}F_{S}~{}r^{2}_{0}\cos\theta~{}d\theta d\phi=0,$ (5)

where $S_{oa}$ is the total ocean surface and $r_{0}$ the radius of the Earth.

  * (iii)

In the results below, we will consider two cases of compensation: (a) global
compensation, i.e. $C$ is the global ocean domain (Fig. 1a) as in Dijkstra
(2007) and (b) $C$ is the Pacific domain (Fig. 1b), so there is no
compensation over the Atlantic. In each case, for different (but fixed) values
of $\gamma_{I}$, a branch of steady solutions versus $\gamma_{A}$ is
calculated under the freshwater forcing (4), starting from the solutions
determined under (i) for $\gamma_{A}=\gamma_{I}=0$.

##  3 Results

In the results below, we concentrate on the bifurcation diagrams and
freshwater and salt balances. Plots of the typical AMOC patterns for slightly
different parameter values can be found in Dijkstra (2007).

###  3.1 Global Compensation

The bifurcation diagram for the global compensation case (Fig. 1a) with
$\gamma_{I}=0$ (no Indian Ocean freshwater flux bias), where the maximum AMOC
strength below 1000 m ($\Psi_{A}$) is plotted versus $\gamma_{A}$ (both in
Sv), is shown as the black curve in Fig. 2a.

Figure 2: _Bifurcation diagrams for the case of global compensation where the
maximum strength of the AMOC below 1000 m ( $\Psi_{A}$) is plotted versus the
strength of the anomalous freshwater forcing $\gamma_{A}$ for different values
of $\gamma_{I}$. Drawn (dashed) curves indicate stable (unstable) branches.
The dots indicate the saddle-node bifurcations. _

With increasing $\gamma_{A}$, stable (black drawn curves) steady states exist
for which the AMOC strength decreases and at $\gamma^{1}_{A}=0.186$ Sv, a
first saddle-node bifurcation $L_{1}$ occurs. With decreasing $\gamma_{A}$, a
branch of unstable steady states (black dashed curves) exists down to a second
saddle-node bifurcation $L_{2}$ at $\gamma^{2}_{A}=0.054$ Sv.

The width of the multi-stable regime, often called the hysteresis width
$\Delta_{H}$, is given by

$\Delta_{H}=|\gamma^{1}_{A}-\gamma^{2}_{A}|.$ (6)

For the case $\gamma_{I}=0$, we find $\Delta_{H}=0.132$ Sv in this model. In
typical quasi-equilibrium model studies Rahmstorf et al. (2005), where
$\gamma_{A}$ is varied with about 0.05 Sv/1000 years, the width is typically
overestimated. With continuation methods, as used here, one is able to
determine the hysteresis width very accurately as the values of
$\gamma^{1,2}_{A}$ are computed explicitly.

With increasing values of $\gamma_{I}$ (adding fresh water over the Indian
Ocean) both saddle node-bifurcations $L_{1}$ and $L_{2}$ move to larger values
of $\gamma_{A}$ (Fig. 2) indicating that the multi-stable regime occurs for
higher Atlantic freshwater forcing. Hence, with the Indian Ocean freshwater
flux bias, it is less likely that GCM having a ‘Levitus-like’ surface salinity
are in a multi-stable regime (than models without such a bias). The width of
the multi-stable regime versus $\gamma_{I}$ does not change much for these
$\gamma_{I}$ values; $\Delta_{H}=0.137$ Sv and $\Delta_{H}=0.138$ Sv are found
for $\gamma_{I}=0.37$ Sv and $\gamma_{I}=0.52$ Sv, respectively.

To understand the shift of the branches in the bifurcation diagram, we
consider the Atlantic freshwater transport by the AMOC and the gyres.
Following De Vries and Weber (2005), the quantities $F_{ov}$ (the overturning
component) and $F_{az}$ (the azonal component) are computed as

$F_{ov}(\theta)&=-\frac{r_{0}}{S_{0}}\int_{S_{\theta}}\bar{v}(<S>-S_{0})~{}dz~{};\\\
F_{az}(\theta)=-\frac{r_{0}}{S_{0}}\int_{S_{\theta}}\overline{v^{\prime}S^{\prime}}~{}dz.$
(7)

where $S_{0}=35$ psu is a reference salinity, $r_{0}$ is the radius of the
Earth, and $S_{\theta}$ is the boundary (longitude, depth) at latitude
$\theta$. Here, the quantities $\bar{v}$, $\bar{S}$, $<v>$ and $<S>$ are given
by

$\bar{v}&=\int
v\cos\theta~{}d\phi~{};~{}<v>=\frac{\bar{v}}{\int\cos\theta~{}d\phi},\\\
\bar{S}=\int
S\cos\theta~{}d\phi~{};~{}<S>=\frac{\bar{S}}{\int\cos\theta~{}d\phi},$ (8)

and $v^{\prime}=v-<v>$ and $S^{\prime}=S-<S>$. The physical meaning of these
quantities is extensively discussed in De Vries and Weber (2005) and Dijkstra
(2007).

The existence of the saddle-node bifurcations $L_{1}$ and $L_{2}$ can be
connected to the behavior of $F_{ovS}$. In simple box models Cessi (1994);
Rahmstorf (1996), the saddle-node bifurcation $L_{1}$ at $\gamma^{1}_{A}$ is
related to a minimum in $F_{ovS}$. In our model, this is more complicated as
there is also a gyre-driven freshwater transport, and the $F_{ovS}$ minimum is
only approximate. The saddle-node bifurcation $L_{2}$ at $\gamma^{2}_{A}$ near
to a zero of $F_{ovS}$ along the upper branch of the AMOC, as discussed at
length in Dijkstra (2007). So to explain the shift in positions in the saddle-
node bifurcations we focus on the behaviour of $F_{ovS}$, $F_{ovN}$ while also
monitoring their difference $\Delta F_{ov}=F_{ovS}-F_{ovN}$ and the associated
behaviour of the freshwater transports by the gyres.

The results for $F_{ovS}=F_{ov}(35^{\circ}$S), $F_{ovN}=F_{ov}(60^{\circ}$N)
and $\Delta F_{ov}$ are shown in Fig. 3a, with the case $\gamma_{I}=0$
($\gamma_{I}=0.37$ Sv) as drawn (dashed) curves.

(a)

(b)

Figure 3: _Global compensation case. (a) Values of the AMOC induced freshwater
transport at the southern boundary (35 ∘S) of the Atlantic ($F_{ovS}$), the
northern boundary (60∘N) of the Atlantic ($F_{ovN}$), and their difference
$\Delta F_{ov}$. (b) Same as (a) but for the azonal transport ($F_{azS}$,
$F_{azN}$ and $\Delta F_{az}$). The drawn curves are for $\gamma_{I}=0$,
corresponding to the black curve in Fig. 2. The dashed curves are for
$\gamma_{I}=0.37$ Sv, corresponding to the blue curve in Fig. 2. _

For the chosen northern latitude (60∘N), the freshwater flux $F_{ovN}$ is
negative and the AMOC transports freshwater southwards for all values of
$\gamma_{A}$. For $\gamma_{A}=0$, the AMOC transports freshwater northwards at
35∘S as $F_{ovS}>0$. With increasing $\gamma_{A}$, $F_{ovS}$ decreases and,
for $\gamma_{I}=0$, becomes negative close to the saddle-node bifurcation
$L_{2}$ at $\gamma^{2}_{A}$ (indicated by the thin drawn vertical line). For
$\gamma_{I}=0.37$ Sv, the value of $F_{ovS}$ at $\gamma_{A}=0$ is much larger
(than that for $\gamma_{I}=0$). Because $F_{ovS}$ decreases with $\gamma_{A}$
at the same rate as for $\gamma_{I}=0$, the location where $F_{ovS}\approx 0$
and hence the bifurcation $\gamma^{2}_{A}$ occurs at a larger value of
$\gamma_{A}$. Also the location where $F_{ovS}$ obtains its minimum, and hence
the position of $L_{1}$ (at $\gamma^{1}_{A}$) shifts to the right. The gyre
transport changes $F_{azS}$, $F_{azN}$ and $\Delta F_{az}$ with $\gamma_{A}$
are shown in Fig. 3b and do not change much with $\gamma_{I}$.

As shown in Dijkstra (2007), the fully-implicit model allows for a closed salt
balance over the Atlantic from $\theta_{s}=35^{\circ}$S to
$\theta_{n}=60^{\circ}$N from which changes in advective, diffusive and
surface contributions can be determined. The terms in this balance are shown
in Fig. 4 for both cases $\gamma_{I}=0$ and $\gamma_{I}=0.37$ Sv. Expressions
for these terms were presented in Dijkstra (2007), but are repeated here for
convenience, i.e.

$\Phi^{s}&=\int_{S_{oa}}S_{0}F_{S}~{}r^{2}_{0}\cos\theta~{}d\phi d\theta,\\\
\Phi^{a}(\theta)=-\int_{S_{\theta}}vS~{}r_{0}\cos\theta~{}d\phi dz,\\\
\Phi^{d}(\theta)=\int_{S_{\theta}}K_{H}\frac{\partial
S}{\partial\theta}~{}\cos\theta~{}d\phi dz,$ (9)

where $\Phi^{a}$ and $\Phi^{d}$ are the advective and diffusive fluxes through
the boundary $S_{\theta}$, respectively. The overall balance is given by

$\Phi^{b}=\Phi^{a}(\theta_{n})-\Phi^{d}(\theta_{n})-\Phi^{a}(\theta_{s})+\Phi^{d}(\theta_{s})-\Phi^{s}$
(10)

Indeed, the term $\Phi^{b}$ is much smaller than the individual terms (black
curves in Fig. 4) giving a nearly closed salt balance over the Atlantic basin
for all values of the parameters.

First consider the case $\gamma_{I}=0$ (drawn curves in Fig. 4) and the upper
branch in the bifurcation diagram up to $L_{1}$. For $\gamma_{A}=0$, the
surface (virtual) salt flux is approximately balanced by the fluxes at the
southern boundary. The surface evaporation is larger than the precipitation
($\Phi^{s}>0$) and this salt is transported out of the Atlantic basin at the
southern boundary ($\Phi^{a}(\theta_{s})<0$). The fact that the diffusive flux
is relatively large here (compared to typical GCMs) is that the model has a
coarse resolution so needs a relatively high horizontal diffusivity to prevent
wiggles to occur near the boundaries. The salt fluxes at the northern boundary
are less important and the respective components are about a factor 2 to 4
smaller than at the southern boundary.

With increasing $\gamma_{A}$, the surface salt flux decreases as freshwater is
put into the North Atlantic. The diffusive salt transports do not respond but
the southward advective salt transport at $\theta_{s}$ weakens and eventually
changes sign near $\gamma_{A}=0.1$ Sv. As the diffusive flux is directed to
transport salt into the basin at the southern boundary, the value of
$\gamma_{A}$ where the sign change in salt transport occurs is around $0.06$
Sv (Fig. 4). Note that the gyre and AMOC components cannot be distinguished in
the advective fluxes.

When $\gamma_{I}=0.37$ Sv (dashed curves in Fig. 4), the surface salt flux
increases for $\gamma_{A}=0$ compared to the case $\gamma_{I}=0$. This is due
to the global compensation as a negative salt flux in the Indian Ocean is
compensated by a positive one over part of the Atlantic. Hence, the curve for
$\Phi^{s}$ shifts upwards and so the compensating advective flux at the
southern boundary shifts downwards. A second effect is that the changed
surface freshwater flux pattern leads to a modified salinity distribution in
the Atlantic. This increases the AMOC (Fig. 2) strength and hence also its
salt transport out of the basin at the southern boundary.

Figure 4: _Global compensation case. Terms in the integrated salt balance
(given in Sv) over the Atlantic basin over the upper branches in Fig. 2 (up to
$\gamma^{1}_{A}$), with expressions of the terms as indicated in (9) and (10)
for $\theta_{n}=60^{\circ}$N and $\theta_{s}=35^{\circ}$S. The drawn curves
are for the case $\gamma_{I}=0$ and the dashed ones for $\gamma_{I}=0.37$ Sv.
_

Because the diffusive fluxes are not much affected by the Indian Ocean
freshwater input, the fluxes $\Phi^{s}$ and $\Phi^{a}(\theta_{s})$ change with
$\gamma_{A}$ in the same way as for the case $\gamma_{I}=0$. The starting
value of $\Phi^{s}$ at $\gamma_{A}=0$ is now larger and it takes a larger
value of $\gamma_{A}$ to change the sign of the freshwater flux at the
southern boundary and to reach a minimum in this quantity. Hence the saddle-
node bifurcations $L_{1}$ and $L_{2}$ shifts to larger values of $\gamma_{A}$.

###  3.2 Pacific Compensation

Since the global compensation has a substantial influence on the position of
the saddle-node bifurcations (Fig. 2), we now consider the case where
compensation is only over the Pacific domain as indicated in Fig. 1b. The
bifurcation diagrams in this case are, for different values of $\gamma_{I}$,
shown in Fig. 5. The shift of the saddle-nodes to larger values of
$\gamma_{A}$ is much smaller than for the global compensation case (Fig. 2).
The hysteresis width itself for $\gamma_{I}=0$ with $\Delta_{H}=0.103$ Sv for
the Pacific compensation case is a bit smaller than for the global
compensation case ($\Delta_{H}=0.138$ Sv). This width is only slightly larger
for the case $\gamma_{I}=0.52$ Sv, i.e. $\Delta_{H}=0.111$ Sv.

Figure 5: _Bifurcation diagram for the case of Pacific compensation where the
maximum strength of the AMOC below 1000 m ( $\Psi_{A}$) is plotted versus the
strength of the anomalous Atlantic freshwater forcing $\gamma_{A}$ for
different values of $\gamma_{I}$. The dots indicate the saddle-node
bifurcations. _

For the analysis of the freshwater and salt balances, we choose the larger
value ($\gamma_{I}=0.52$ Sv) instead of $\gamma_{I}=0.37$ Sv (used in the
global compensation case), as the differences are more clearly visible. The
freshwater transports by the AMOC and by the gyres (Fig. 6) show that for
$\gamma_{A}=0$, $F_{ov}$ is larger for $\gamma_{I}=0.52$ Sv (Fig. 6a),
compared to the case $\gamma_{I}=0.0$ Sv. Hence also changes in the freshwater
balance are induced in the Atlantic, but both the direct effect of
compensation of the surface freshwater flux and the secondary effect of an
AMOC increase are much smaller. Of course, in such a diffusive and viscous
model, the Agulhas retroflection is in a diffusive retroflection regime
Dijkstra and De Ruijter (2001) but there is additional freshwater transport
from the Indian to the Atlantic when $\gamma_{I}>0$. As the wind-driven
freshwater transport does not change much with $\gamma_{A}$ (Fig. 6b), the
AMOC transports more freshwater into the basin and hence $F_{ovS}$ becomes
more positive.

(a)

(b)

Figure 6: _Pacific compensation case. (a) Values of the AMOC induced
freshwater transport at the southern boundary (35 ∘S) of the Atlantic
($F_{ovS}$), the northern boundary (60∘N) of the Atlantic ($F_{ovN}$), and
their difference ($\Delta F_{ov}$). (b) Same as (a) but for the azonal
transport ($F_{azS}$, $F_{azN}$ and $\Delta F_{az}$). The drawn curves are for
$\gamma_{I}=0$, corresponding to the black curve in Fig. 5. The dashed curves
are for $\gamma_{I}=0.52$ Sv, corresponding to the red curves in Fig. 5. _

In the Atlantic salt balance (Fig. 7) the surface salt flux is indeed the same
for both values of $\gamma_{I}$, because there is no compensation anymore over
the Atlantic. The advective transport of salt becomes slightly more negative
over the southern boundary for $\gamma_{I}=0.52$ Sv indicating that indeed a
small amount of freshwater is transported into the Atlantic basin. Also the
diffusive transport of salt (into the basin) decreases and this approximately
balances the advective contribution. As the advective contribution is much
smaller than the compensation contribution in the global compensation case,
the shift of the saddle-node bifurcations to larger values of $\gamma_{A}$ is
much smaller.

Figure 7: _Pacific compensation case. Terms in the integrated salt balance
over the Atlantic basin over the upper branch in Fig. 5, with expressions of
the terms as indicated in the main text. The drawn curves are for the case
$\gamma_{I}=0$ and the dashed ones for $\gamma_{I}=0.52$ Sv. _

##  4 Summary and Discussion

Following earlier studies on CMIP3 and CMIP5 models Drijfhout et al. (2013);
Mecking et al. (2016), also many CMIP6 models have large biases in surface
freshwater fluxes which lead to an AMOC with an Atlantic freshwater transport
that is in disagreement with observations van Westen and Dijkstra (2023). The
most important model bias is in the Atlantic Surface Water properties, which
arises from deficiencies in the surface freshwater flux over the Indian Ocean.
A second bias is in the properties in the North Atlantic Deep Water and arises
through deficiencies in the freshwater flux over the Atlantic Subpolar Gyre
region van Westen and Dijkstra (2023).

In this paper, we have addressed to effects of the freshwater flux bias in the
Indian Ocean on the multiple equilibrium regime of the AMOC using the fully
implicit global ocean-atmosphere model of Dijkstra and Weijer (2005), for
which explicit bifurcation diagrams can be computed. This is a fantastic
capability as both stable and unstable steady states can be computed, and the
width of the multiple equilibrium regime can be determined accurately.
However, this can only be done with quite a simplified global model, with
relatively low resolution and hence being more viscous and diffusive than
current, even low-resolution, ocean models. The atmosphere model is only an
energy balance model with a prescribed freshwater forcing. In terms of
quantitative results on changes of the bifurcation diagrams, the model is
probably not that useful.

Qualitatively, however, the model provides very useful information, as it
indicates that the multiple equilibrium regime does not disappear due to the
Indian Ocean freshwater flux biases, but that it shifts to higher values of
the North Atlantic anomalous freshwater flux. This shift is dependent on the
way the latter flux is compensated. There are two mechanisms which are
responsible for these shifts: surface salinity patterns in the Atlantic depend
on the compensation of the Indian Ocean bias and, in both cases considered
here, lead to a slight increase in the AMOC. In the Levitus background state
(with $F_{ovS}>0$), this leads to a larger transport of salt out of the
Atlantic basin. Hence, a larger anomalous North Atlantic surface freshwater
flux is needed to activate the salt advection feedback, i.e., to a situation
where the AMOC exports fresh water. When there is compensation of the Indian
Ocean freshwater flux in the Atlantic, the Atlantic becomes saltier and hence
the AMOC transports more salt out of the basin. This leads to an additional,
and here larger, shift of the bifurcation diagram compared to when there is no
compensation in the Atlantic.

The mechanisms identified are useful for interpreting results from GCMs and
for designing new simulations with these models. First it shows that biases in
freshwater flux lead to shifts in the bifurcation diagram which would imply
that such multiple equilibrium regimes (and hence AMOC collapses) could exist
in these models but are located in a parameter regime where one would not
normally perform simulations. Second, it provides a hint why efforts to find
AMOC collapses in these models may not have been successful Mecking et al.
(2016); Jackson and Wood (2018b, a). Even an enormous freshwater input in a
parameter regime which is outside the multiple equilibrium regime would only
lead to a weakened (but no collapsed) AMOC.

This also suggests several ways forward to find an AMOC collapse in state-of-
the-art models. One either performs a long quasi-equilibrium simulation up to
very large freshwater flux input such as in the FAMOUS model Hawkins et al.
(2011) to find the collapse (this may take a few thousand years of simulation,
so is expensive). In doing this, it is better to compensate outside of the
Atlantic when using surface fluxes, as the latter will introduce an additional
shift and so one has to integrate longer to find an AMOC collapse. Such
compensation procedures (including compensation over the volume) are now also
more common in GCMs Jackson et al. (2022). The alternative is to address and
improve the biases in the atmospheric components of the models, but this is
not an easy issue. The origin of these biases may even be a coupled problem,
as the bias strength is positively correlated with the AMOC strength van
Westen and Dijkstra (2023).

We hope that this study will motivate the design of new AMOC hosing
simulations to find AMOC collapses in state-of-the-art GCMs. Detection of such
a collapse would have a large impact on climate change research and
probability estimates of AMOC tipping under global warming would likely need
to be revised.

Acknowledgements

H.A.D. and R.M.v.W. are funded by the European Research Council through the
ERC-AdG project TAOC (project 101055096). The authors thank Dr. Fred Wubs
(University of Groningen, NL) for helping with revising the old Fortran code
of the model Dijkstra (2007) used in this paper.

## References

  * Armstrong McKay et al. (2022) Armstrong McKay, D. I., Staal, A., Abrams, J. F., Winkelmann, R., Sakschewski, B., Loriani, S., Fetzer, I., Cornell, S. E., Rockström, J., and Lenton, T. M. (2022).  Exceeding 1.5 C global warming could trigger multiple climate tipping points.  Science, 377(6611), eabn7950.
  * Bryan and Lewis (1979) Bryan, K. and Lewis, L. J. (1979).  A water mass model of the world ocean.  J. Geophys. Res., 84, 2503–2517.
  * Bryden et al. (2011) Bryden, H. L., King, B. A., and McCarthy, G. D. (2011).  South Atlantic overturning circulation at 24S.  Journal of Marine Research, 69, 38–55.
  * Cessi (1994) Cessi, P. (1994).  A simple box model of stochastically forced thermohaline flow.  J. Phys. Oceanogr., 24, 1911–1920.
  * Cimatoribus et al. (2012) Cimatoribus, A., Drijfhout, S., Toom, M., and Dijkstra, H. (2012).  Sensitivity of the Atlantic meridional overturning circulation to South Atlantic freshwater anomalies.  Climate Dynamics, 39, 2291 – 2306.
  * De Niet et al. (2007) De Niet, A., Wubs, F., van Scheltinga, A. T., and Dijkstra, H. A. (2007).  A tailored solver for bifurcation analysis of ocean-climate models.  Journal of Computational Physics, 227(1), 654–679.
  * De Vries and Weber (2005) De Vries, P. and Weber, S. L. (2005).  The Atlantic freshwater budget as a diagnostic for the existence of a stable shut down of the meridional overturning circulation.  Geophys. Res. Letters, 32, L09606.
  * Dijkstra (2005) Dijkstra, H. A. (2005).  Nonlinear Physical Oceanography: A Dynamical Systems Approach to the Large Scale Ocean Circulation and El Niño, 2nd Revised and Enlarged edition. Springer, New York, 532 pp.
  * Dijkstra (2007) Dijkstra, H. A. (2007).  Characterization of the multiple equilibria regime in a global ocean model.  Tellus, 59A, 695–705.
  * Dijkstra and De Ruijter (2001) Dijkstra, H. A. and De Ruijter, W. P. M. (2001).  On the physics of the Agulhas current: Steady retroflection regimes.  J. Phys. Oceanogr., 31, 2971–2985.
  * Dijkstra and Weijer (2005) Dijkstra, H. A. and Weijer, W. (2005).  Stability of the global ocean circulation: basic bifurcation diagrams.  J. Phys. Oceanogr., 35, 933–948.
  * Drijfhout et al. (2013) Drijfhout, S., Gleeson, E., Dijkstra, H. A., and Levina, V. (2013).  Spontaneous abrupt climate change due to an atmospheric blocking– sea-ice–ocean feedback in an unforced climate model simulation.  Proceedings National Acad. Sciences.
  * England (1993) England, M. H. (1993).  Representing the global-scale water masses in ocean general circulations models.  J. Phys. Oceanogr., 23, 1523–1552.
  * Gent (2018) Gent, P. R. (2018).  A commentary on the atlantic meridional overturning circulation stability in climate models.  Ocean Modelling, 122, 57–66.
  * Hawkins et al. (2011) Hawkins, E., Smith, R. S., Allison, L. C., Gregory, J. M., Woollings, T. J., Pohlmann, H., and De Cuevas, B. (2011).  Bistability of the Atlantic overturning circulation in a global climate model and links to ocean freshwater transport.  Geophysical Research Letters, 38(10), L10605.
  * Hu et al. (2012) Hu, A., Meehl, G. A., Han, W., Timmermann, A., Otto-Bliesner, B., Liu, Z., Washington, W. M., Large, W., Abe-Ouchi, A., Kimoto, M., Lambeck, K., and Wu, B. (2012).  Role of the Bering Strait on the hysteresis of the ocean conveyor belt circulation and glacial climate stability.  Proceedings of the National Academy of Sciences, 109(17), 6417–6422.
  * Huisman et al. (2010) Huisman, S. E., den Toom, M., Dijkstra, H. A., and Drijfhout, S. (2010).  An Indicator of the Multiple Equilibria Regime of the Atlantic Meridional Overturning Circulation.  Journal Of Physical Oceanography, 40(3), 551–567.
  * Jackson and Wood (2018a) Jackson, L. C. and Wood, R. A. (2018a).  Hysteresis and Resilience of the AMOC in an Eddy‐Permitting GCM.  Geophysical Research Letters, 45(16), 8547 – 8556.
  * Jackson and Wood (2018b) Jackson, L. C. and Wood, R. A. (2018b).  Timescales of AMOC decline in response to fresh water forcing.  Climate Dynamics, 51(4), 1333 – 1350.
  * Jackson et al. (2015) Jackson, L. C., Kahana, R., Graham, T., Ringer, M. A., Woollings, T., Mecking, J. V., and Wood, R. A. (2015).  Global and European climate impacts of a slowdown of the AMOC in a high resolution GCM.  Climate Dynamics, 45(11), 3299 – 3316.
  * Jackson et al. (2022) Jackson, L. C., Asenjo, E. A. d., Bellomo, K., Danabasoglu, G., Haak, H., Hu, A., Jungclaus, J., Lee, W., Meccia, V. L., Saenko, O., Shao, A., and Swingedouw, D. (2022).  Understanding AMOC stability: the North Atlantic Hosing Model Intercomparison Project.  Geoscientific Model Development Discussions, 2022, 1–32.
  * Keller (1977) Keller, H. B. (1977).  Numerical solution of bifurcation and nonlinear eigenvalue problems.  In P. H. Rabinowitz, editor, Applications of Bifurcation Theory. Academic Press, New York, U.S.A.
  * Lenton et al. (2008) Lenton, T. M., Held, H., Kriegler, E., Hall, J. W., Lucht, W., Rahmstorf, S., and Schellnhuber, H. J. (2008).  Tipping elements in the Earth’s climate system.  Proceedings of the National Academy of Sciences of the United States of America, 105(6), 1786–93.
  * Levitus (1994) Levitus, S. (1994).  World Ocean Atlas 1994, Volume 4: Temperature.  NOAA/NESDIS E, US Department of Commerce, Washington DC, OC21, 1–117.
  * Liu et al. (2017) Liu, W., Xie, S.-P., Liu, Z., and Zhu, J. (2017).  Overlooked possibility of a collapsed Atlantic Meridional Overturning Circulation in warming climate.  Science Advances, 3(1), e1601666.
  * Liu et al. (2020) Liu, W., Fedorov, A. V., Xie, S.-P., and Hu, S. (2020).  Climate impacts of a weakened Atlantic Meridional Overturning Circulation in a warming climate.  Science Advances, 6(26), eaaz4876.
  * Lohmann et al. (2021) Lohmann, J., Castellana, D., Ditlevsen, P. D., and Dijkstra, H. A. (2021).  Abrupt climate change as rate-dependent cascading tipping point.  Earth System Dynamics Discussions, 2021, 1–25.
  * Marotzke (2000) Marotzke, J. (2000).  Abrupt climate change and thermohaline circulation: Mechanisms and predictability.  Proc. Natl. Acad. Sci., 97, 1347–1350.
  * Mecking et al. (2016) Mecking, J. V., Drijfhout, S. S., Jackson, L. C., and Graham, T. (2016).  Stable AMOC off state in an eddy-permitting coupled climate model.  Climate Dynamics, 47(7), 2455 – 2470.
  * Mulder et al. (2021) Mulder, T. E., Goelzer, H., Wubs, F. W., and Dijkstra, H. A. (2021).  Snowball Earth Bifurcations in a Fully-Implicit Earth System Model.  International Journal of Bifurcation and Chaos, 31(06), 2130017\.
  * Orihuela-Pinto et al. (2022) Orihuela-Pinto, B., England, M. H., and Taschetto, A. S. (2022).  Interbasin and interhemispheric impacts of a collapsed atlantic overturning circulation.  Nature Climate Change, 12(6), 558–565.
  * Rahmstorf (1996) Rahmstorf, S. (1996).  On the freshwater forcing and transport of the Atlantic thermohaline circulation.  Climate Dynamics, 12(12), 799–811.
  * Rahmstorf et al. (2005) Rahmstorf, S., Crucifix, M., Ganopolski, A., Goosse, H., Kamenkovich, I., Knutti, R., Lohmann, G., March, R., Mysak, L., Wang, Z., and Weaver, A. J. (2005).  Thermohaline circulation hysteresis: a model intercomparison.  Geophys. Res. Letters, L23605, doi:0.1029/2005GLO23655, 1–5.
  * Stommel (1961) Stommel, H. (1961).  Thermohaline convection with two stable regimes of flow.  Tellus, 2, 244–230.
  * Stouffer et al. (2006) Stouffer, R. J., Yin, J., Gregory, J. M., Dixon, K. W., Spelman, M. J., Hurlin, W., Weaver, A. J., Eby, M., Flato, G. M., Hasumi, H., Hu, A., Jungclaus, J. H., Kamenkovich, I. V., Levermann, A., Montoya, M., Murakami, S., Nawrath, S., Oka, A., Peltier, W. R., Robitaille, D. Y., Sokolov, A. P., Vettoretti, G., and Weber, S. L. (2006).  Investigating the causes of the response of the thermohaline circulation to past and future climate changes.  Journal of Climate, 19, 1365–1387.
  * Toom et al. (2012) Toom, M. D., Dijkstra, H. a., Cimatoribus, A. a., and Drijfhout, S. S. (2012).  Effect of Atmospheric Feedbacks on the Stability of the Atlantic Meridional Overturning Circulation.  Journal of Climate, 25(12), 4081–4096.
  * Trenberth et al. (1989) Trenberth, K. E., Olson, J. G., and Large, W. G. (1989).  A global ocean wind stress climatology based on ECMWF analyses.  Technical report, National Center for Atmospheric Research, Boulder, CO, U.S.A.
  * van Westen and Dijkstra (2023) van Westen, R. M. and Dijkstra, H. A. (2023).  Persistent Climate Model Biases in the Atlantic Ocean’s Freshwater Transport.  Ocean Science Discussions.  submitted.
  * Vellinga et al. (2002) Vellinga, M., Wood, R. A., and Gregory, J. M. (2002).  Processes governing the recovery of a perturbed thermohaline circulation in HadCM3.  J. Climate, 15, 764–780.
  * Weijer et al. (2019) Weijer, W., Cheng, W., Drijfhout, S. S., Fedorov, A. V., Hu, A., Jackson, L. C., Liu, W., McDonagh, E. L., Mecking, J. V., and Zhang, J. (2019).  Stability of the Atlantic Meridional Overturning Circulation: A Review and Synthesis.  Journal Of Geophysical Research-Oceans, 124(8), 5336 – 5375\.


# Digital signature schemes using non-square matrices or scrap automorphisms

Jiale Chen  Department of Mathematics, The City College of New York, New York,
NY 10031<EMAIL_ADDRESS>,  Dima Grigoriev  CNRS, Mathématiques,
Université de Lille, 59655, Villeneuve d’Ascq, France  Dmitry.Grigoryev@univ-
lille.fr and  Vladimir Shpilrain  Department of Mathematics, The City College
of New York, New York, NY 10031<EMAIL_ADDRESS>
###### Abstract.

We offer two very transparent digital signature schemes: one using non-square
matrices and the other using scrap automorphisms. The former can be easily
converted to a public key encryption scheme.

##  1\. Introduction

Due to the concern that if large-scale quantum computers are ever built, they
will compromise the security of many commonly used cryptographic algorithms,
NIST had begun in 2016 a process to develop new cryptography standards and, in
particular, solicited proposals for new digital signature schemes [7]
resistant to attacks by known quantum algorithms, such as e.g. [9].

One possible way to avoid quantum attacks based on solving the hidden subgroup
problem (including the attacks in [9]) is not to use one-way functions on sets
that have an obvious (semi)group structure. In this paper, we offer two
digital signature schemes that are similar in spirit. What they have in common
is the following idea. There is a function $F:X\to Y$ that is not necessarily
one-way. However, when only a part (call it $\bar{Y}$) of $Y$ is published,
the restriction $\bar{F}:\bar{X}\to\bar{Y}$ becomes much harder to invert and
therefore becomes a candidate for a one-way function.

One of the simplest instantiations of this general idea is employed in our
first digital signature scheme. There, $X$ is an invertible square matrix over
a commutative ring and $Y$ is $X^{-1}$. The function $F$ therefore takes a
square matrix as the input and outputs the inverse matrix. This is often not a
one-way function. However, if we only publish a subset of columns of $X^{-1}$,
then finding the corresponding subset of rows of $X$ appears to be a hard
problem. The corresponding digital signature scheme is described in Section 2.
We note that the verification part in this scheme is very efficient since it
amounts to multiplying a vector of sparse polynomials by a matrix whose
entries are sparse polynomials (over $\mathbb{Z}_{q}$, for a small $q$).

For our second scheme, we consider an automorphism $\alpha:K\to K$ of an
algebra $K$ of polynomials in $n$ variables $x_{1},\ldots,x_{n}$ over a ring
$R$. Such an automorphism can be given by $n$ polynomials $\alpha(x_{i})$. It
is arguable whether or not computing $\alpha^{-1}$ from all $\alpha(x_{i})$ is
a hard problem (this might depend on the ring $R$). Some of the known methods
use Gröbner bases techniques (see e.g. [10]) that are not known to be
efficient. A more efficient algorithm was reported in [2], but it is still
exponential-time in the number of variables, although polynomial-time in the
number of polynomials and their degrees. We note that there was a proposal
some time ago for a public key cryptosystem based on the alleged hardness of
inverting polynomial automorphisms [6], but it was subsequently attacked in
[3].

In the present paper, we suggest to publish not all of the $\alpha(x_{i})$,
but only some of them, with the idea that inverting such a “scrap”
automorphism should be harder than inverting the “whole” $\alpha$. To the best
of our knowledge, the only known way to invert a “scrap” automorphism is the
“brute force” trial-and-error method.

Our digital signature scheme based on “scrap” automorphisms is described in
Section 11.

Finally, we note that our digital signature scheme based on scrap matrices can
be easily converted to a public key encryption scheme, see Section 10. We do
not make any security claims about this scheme though since that would require
a separate security analysis, and this is not in the scope of the present
paper.

##  2\. Non-square matrices: scheme description

Let $K=\mathbb{Z}_{q}[x_{1},\ldots,x_{n}]$ denote the algebra of polynomials
in $n$ variables over the ring $\mathbb{Z}_{q}$ of integers modulo $q$, where
an integer $q$ is not necessarily a prime.

The signature scheme is as follows.

Public:

– $k\times l$ left invertible matrix $M$, with $k>l$, whose entries are sparse
polynomials from the algebra $K$.
– a hash function $H$ and a (deterministic) procedure for converting values of
$H$ to vectors of sparse polynomials from the algebra $K$.

Private: $l\times k$ right invertible matrix $L$ over $K$, such that $LM$ is
the $l\times l$ identity matrix.

Signing a message $m$:

  1. (1)

Apply a hash function $H$ (e.g., SHA-512) to $m$. Convert $H(m)$ to a vector
${\bf U}=(P_{1},\ldots,P_{l})$ of $l$ (sparse) polynomials from the algebra
$K$ using a deterministic public procedure.

  2. (2)

Multiply the vector ${\bf U}$ by the (private) matrix $L$ on the right to get
a vector ${\bf V}={\bf U}L=(Q_{1},\ldots,Q_{k})$ of $k$ polynomials from $K$.

  3. (3)

The signature is the vector ${\bf V}$.

Verification.

1. The verifier computes the hash $H(m)$ and converts $H(m)$ to a vector ${\bf U}=(P_{1},\ldots,P_{l})$ of $l$ (sparse) polynomials using a deterministic public procedure.

2. The verifier multiplies the signature vector ${\bf V}$ by the public matrix $M$ on the right to get a vector ${\bf W}$.

3. The signature is accepted if and only if ${\bf W}={\bf U}$.

Correctness is obvious since ${\bf W}={\bf V}M=({\bf U}L)M={\bf U}(LM)={\bf
U}$.

######  Remark 1.

One can use a faster, numerical, variant of the signature verification.
Namely, instead of computing a vector ${\bf W}$ of polynomials, the verifier
can assign random numerical values (from $\mathbb{Z}_{q}$) to all variables
$x_{i}$ in the vectors ${\bf U}$, ${\bf V}$ and the matrix $M$, and then do
computation in $\mathbb{Z}_{q}$. To make the probability of the “false
positive” decision negligible, assigning random numerical values should be
done several times.

##  3\. Key generation

In this section, we answer crucial questions on key generation in our scheme:
how to generate a (left) invertible $k\times l$ matrix, how to compute its
inverse, and how to generate a random sparse polynomial. We emphasize that in
our scheme, all this is done offline.

###  3.1. Generating a random $t$-sparse polynomial

We call a polynomial $t$-sparse if it has $t$ monomials. To generate a
“random” $t$-sparse polynomial, we first generate $t$ monomials as follows.

1. Bound the degree of a monomial by an integer $b$ (one of the parameters of the scheme).

2. Select, uniformly at random, an integer $d$ between 1 and $b$. This integer will be the degree of our monomial.

3. To select a monomial of degree $d$, we do a selection of $x_{i}$, uniformly at random from $\\{x_{1},\ldots,x_{n}\\}$, $d$ times. Then build the monomial as a product of the selected $x_{i}$.

4. Having selected $t$ monomials like that, build a $t$-sparse polynomial as a linear combination of $t$ selected monomials with random coefficients from $\mathbb{Z}_{q}-\\{0\\}$.

###  3.2. Generating a (left) invertible $k\times l$ matrix

Let $k>l$. To generate a (left) invertible $k\times l$ matrix, one can first
generate an invertible square $k\times k$ matrix and then remove $k-l$
columns, selected at random.

To generate an invertible square $k\times k$ matrix, one can do the following.

1. Generate an upper unitriangular $k\times k$ matrix $U$ with random $t$-sparse polynomials above the diagonal and a lower unitriangular $k\times k$ matrix $L$ with random $t$-sparse polynomials below the diagonal.

2. Make a product of matrices that would have the form $UP_{1}LP_{2}$, where $P_{i}$ are matrices corresponding to random permutations of columns and rows of a $k\times k$ matrix. (Entries of $P_{i}$ are 0s and 1s.)

###  3.3. Computing the (left) inverse of a matrix

Having generated an invertible square matrix $S=UP_{1}LP_{2}$ (Section 3.2),
we compute its inverse as
$S^{-1}=P_{2}^{-1}L^{-1}P_{1}^{-1}U^{-1}$. Computing $P_{i}^{-1}$ is trivial,
and computing the inverse of a unitriangular square matrix $X$ is done by
using the adjugate of $X$. (Note that the determinant of any unitriangular
square matrix is equal to 1.)

Now suppose the (left) invertible $k\times l$ matrix $M$ was obtained from the
square $k\times k$ matrix $S$ by removing $k-l$ columns
$C_{i_{1}},\ldots,C_{i_{k-l}}$. Then, to get a left inverse of $M$, we just
remove the corresponding rows $R_{i_{1}},\ldots,R_{i_{k-l}}$ from $S^{-1}$.

###  3.4. Converting $H(m)$ to a vector of sparse polynomials

We suggest using a hash functions from the SHA-2 family, specifically SHA-512.
We assume the security properties of SHA-512, including collision resistance
and preimage resistance. The following procedure may seem “too ad hoc”, but
this is needed to make sure that the linear space of possible monomials
obtained by this procedure has a large dimension, to prevent the forgery
attack described in Section 7.1.

Let $B=H(m)$ be a bit string of length 512. We will convert $B$ to a vector
${\bf U}=(P_{1},\ldots,P_{l})$ of $l$ (4-sparse) polynomials from the algebra
$K$, where each monomial has degree at most 10. We note that this process is
deterministic. Of course, this is just one of the many possible conversion
procedures; one of our goals here is to only have polynomials $P_{i}$ of low
degree.

(1) Discard several last bits in $B$ to make its length divisible by $l$. With
our suggested parameters (see Section 4), $l=5$, so one can discard, say, the
last 12 bits of $B$.

(2) Split the remaining $500$ bits in two parts: $300$ in one part and $200$
in the other. The part with $300$ bits will determine a map of the set
$\\{x_{1},\ldots,x_{50}\\}$ of variables to the set
$\\{x_{1},\ldots,x_{63}\\}$, as follows.

Each integer from 1 to 50 can be represented in the binary form using 6 bits.
Thus, $300=6\cdot 50$ bits will naturally determine a map $x_{i}\to x_{j}$ for
each $i$ from 1 to 50. For example, if the first 6 bits of these 300 are
011010, this means that $x_{1}$ will be mapped to $x_{26}$ because 011010=26.

This map is not necessarily one-to-one (since some of the 6-bit blocks may be
identical) but this is not a problem. Let $\\{x_{i_{1}},\ldots,x_{i_{50}}\\}$
be the new sequence of variables; some of them may be equal.

(3) Split the remaining $200$ bits in $l=5$ 40-bit blocks. Then split each
40-bit block in 4 parts, 10 bits in each. This will correspond to 4 monomials,
of degree at most 10 each, as follows. Enumerate bits in each 40-bit block by
$x_{i_{1}},\ldots,x_{i_{40}}$ (in this order, going left to right). Now each
block of 10 bits is converted to a monomial that is a product of $x_{i}$
corresponding to the places in the bit string where the bit is “1”. In
particular, each monomial will be of degree at most 10. For example, if the
first 10-bit block is 1001100011, then the corresponding monomial will be
$x_{i_{1}}x_{i_{4}}x_{i_{5}}x_{i_{9}}x_{i_{10}}$.

(4) Combine 4 monomials obtained at Step (3) into a 4-sparse polynomial by
selecting coefficients from $\mathbb{Z}_{q}$ as follows. We are going to use
the string of 12 bits that we have discarded at Step (1). Split it in 4 blocks
of 3 bits. Each block will determine a coefficient at a monomial in a 4-sparse
polynomial by converting a 3-bit binary number to decimal and reducing modulo
6.

##  4\. Suggested parameters

For the hash function $H$, we suggest SHA-512.

For the integer $q$ in $\mathbb{Z}_{q}$, we suggest $q=6$.

For the number $n$ of variables, we suggest $n=64$.

For the dimensions of the matrix $M$, we suggest $k=10,~{}l=5$. However, due
to current software limitations (see Section 6 for more details), our present
implementation uses $k=5,~{}l=3$.

For the number $t$ of monomials in $t$-sparse polynomials that entries of the
unitriangular matrices, we suggest $t=3$.

For the bound $b$ on the degree of monomials in $t$-sparse polynomials, we
suggest $b=3$.

##  5\. What is the hard problem here?

The (computationally) hard problem that we employ in our construction is
finding a left inverse of a given left invertible matrix $M$. That is, solving
the matrix equation $XM=I$, where $X$ is the unknown matrix of given
dimensions.

If matrices in this equation were considered over a field, then this matrix
equation would translate to a system of linear equations (in the entries of
the matrix $X$), and therefore would be easily solvable.

In our situation, where matrices are considered over a polynomial algebra, the
problem can still be reduced to a system of linear equations, this time in the
coefficients of polynomials that are entries of the matrix $X$. If the number
$n$ of variables $x_{i}$ is not too small (we suggest $n=64$), then the number
of different monomials (and therefore the number of unknown coefficients of
polynomials) is quite large (see Section 7), so that the relevant system of
linear equations becomes intractable.

###  5.1. Completing a left invertible matrix to an invertible square matrix

Another possible way to find a left inverse of a given left invertible matrix
$M$ is to complete $M$ to an invertible square matrix by adding more columns,
and then find the inverse of this square matrix (which is relatively easy
since one can use determinants).

Perhaps surprisingly, it was a major open problem (Serre’s problem) in algebra
whether or not any left invertible matrix over a polynomial algebra can be
completed to an invertible square matrix by adding more rows (or columns).
This problem was answered in the affirmative independently by Quillen [8] and
Suslin [11].

The question that matters to us though is that of the computational complexity
of completing $M$ to an invertible square matrix. It was shown in [1] (for
polynomial algebras over an infinite field) and then in [5] (for polynomial
algebras over an arbitrary Noetherian ring) that there is a relevant algorithm
whose complexity is exponential in the square of the number of variables
$x_{i}$. This is yet another reason why the number of variables should not be
too small.

##  6\. Performance and signature size

Disclaimer. We have run into a problem with software for symbolic computation
with multivariate polynomials, especially with matrices over multivariate
polynomials. None of the symbolic computation packages either in Julia or in
ANSI C (Flint) is capable of doing expand-and-simplify if the number of
monomials becomes as large as several hundred. The only software that is
capable of doing that is Maple, but since it is not common to use Maple for
cryptographic purposes, we have decided not to use it; instead, we have
reduced the values of the most impactful parameters, the dimensions of the
public matrix $M$, from the recommended $10\times 5$ to $5\times 3$. With
Flint evolving, it should be soon possible to handle $10\times 5$ matrices
over multivariate polynomials.

For our computer simulations, we used Apple MacBook Pro, M1 CPU (8 Cores), 16
GB RAM computer.

With the suggested parameters, signature verification takes about 0.2 sec on
average, which is not bad, but the signature is rather large, about 4,200
bytes on average.

The size of the private key (the matrix $L$) is about 2,000 bytes, and so is
the size of the public key (the matrix $M$).

We note that we have measured the size of a signature, as well as the size of
private/public keys, as follows. We have counted the total number of variables
that occurred in relevant polynomial(s) and multiplied that number by 7, the
number of bits sufficient to describe the index of any variable (except
$x_{64}$). To that, we added the number of monomials times 2 (the average
number of bits needed to describe a coefficient at a monomial in our
construction(s)).

Performance metrics for given parameter values
---
# monomials in sparse polynomials (parameter $t$) |  max degree of monomials (parameter $b$) |  dimensions of private matrix $L$ |  verification time (sec) |  signature size (Kbytes) |  public key size (Kbytes) |  private key size (Kbytes) |  memory usage (Mbytes)
3 |  3 |  3x5 |  0.2 |  4.2 |  2 |  1.9 |  1.4

##  7\. Linear algebra attacks

One can attempt to recover the private left inverse of the public matrix $M$
by using linear aOne can attempt to recover the private left inverse of the
public matrix $M$ by using linear algebra, as follows. Finding the left
inverse matrix of $M$ is equivalent to solving the matrix equation $XM=I$,
where $X$ is the unknown matrix of given dimensions.

If matrices in this equation were considered over a field, then this matrix
equation would translate to a system of linear equations (in the entries of
the matrix $X$), and therefore would be easily solvable.

In our situation, where matrices are considered over a polynomial algebra, the
problem can still be reduced to a system of linear equations, this time in the
coefficients of polynomials that are entries of the matrix $X$. If the number
$n$ of variables $x_{i}$ is not too small (we suggest $n=64$), then the number
of different monomials (and therefore the number of unknown coefficients of
polynomials) is very large, so that the relevant system of linear equations
becomes intractable.

Note that, although monomials in the public matrix $M$ have degrees bounded by
3, in the left inverse of $M$ there can be monomials of significantly higher
degree. Let us estimate the number of monomials of degree, say, $15$ in 64
variables.

By the formula for the number of combinations with repetitions, the number is
${78\choose 15}\approx 2^{52}.$

Thus, the total number of potentially possible monomials in a $3\times 6$
matrix with sparse polynomial entries is more than $2^{52}\cdot 15>2^{55}$,
and this would be the number of unknowns in a system of linear equations. The
number of equations would be at least as large.

At least $(2^{55})^{2.3}\approx 2^{126}$ arithmetic operations are needed to
solve such a system of linear equations, according to our understanding of the
state-of-the-art in solving systems of linear equations.

###  7.1. Forgery after collecting multiple signatures

A smarter linear algebra attack is targeted at forging signatures without
recovering the private key and is based on collecting multiple vectors ${\bf
U}=(P_{1},\ldots,P_{l})$ and finding a basis of the linear space generated by
these vectors.

Since the transformation ${\bf U}\to{\bf V}={\bf U}L$ is linear, after a basis
of the linear space of vectors ${\bf U}$ is constructed, the forger can
express each newly obtained ${\bf U}$ as a linear combination of previously
collected ones, and then the signature corresponding to this newly obtained
${\bf U}$ will be the same linear combination of the relevant previously
published relevant signatures ${\bf V}$.

To make this attack infeasible, the dimension of the linear space of possible
vectors ${\bf U}$ should be quite large. This dimension is equal to the number
of monomials that can possibly appear in polynomials $P_{1},\ldots,P_{l}$
under the procedure of converting $H(m)$ to a vector of sparse polynomials,
see Section 3.4. The number of monomials of degree $k\leq 64$ in 64 variables
is ${64+k-1\choose k}$, so in our situation, where the degree can be any
number between 1 and 15, the number of possible monomials is at least
${78\choose 15}\approx 2^{52}$. This is about how many linearly independent
signatures a forger has to collect. When collecting, say, 1000 linearly
independent signatures a second, this would take over 50000 years.

##  8\. Security claims

With suggested parameters, the linear algebra “brute force” attack requires at
least $2^{126}$ arithmetic operations in $\mathbb{Z}_{6}$ to solve the
relevant system of linear equations.

There could be ad hoc attacks on the public key aiming at recovering some of
the entries of the private matrix, but recovering only some of the entries
cannot make the probability of passing verification non-negligible because
getting even just a couple of entries (out of 15) wrong will not increase
chances of passing verification.

Forgery without getting a hold of the private key requires accumulation of at
least $2^{52}$ linearly independent signatures.

As for quantum security, it appears that no known quantum algorithms are
applicable in our situation since there are no abelian (semi)groups in play
here that could be exploited by one of the known quantum attacks.

##  9\. Advantages and limitations of the matrix digital signature scheme

###  9.1. Advantages

1. Unparalleled conceptual simplicity.

2. Simplicity of design. Signing a message amounts to just a matrix-vector multiplication, and so does verifying a signature.

3. Bringing into play an interesting hard problem – inverting a non-square matrix over multivariate polynomials. This problem has not been considered in the cryptographic context before, to the best of our knowledge.

4. Efficiency of the signature verification (about 0.2 sec on average with parameters in the current implementation). To be fair, if we did not have to downsize dimensions of the public matrix from $10\times 5$ to $5\times 3$ (see the Disclaimer in Section 6), then we project the signature verification time would be close to 1 sec.

###  9.2. Limitations

1. Linearity of the signing function. To avoid linear algebra attacks, one has to use a large number of variables in the polynomial algebra, and this can affect the signature size (but not so much the verification time).

2. Large signature size.

##  10\. Public key encryption

Our first digital signature scheme (Section 2) can be easily converted to a
public key encryption scheme, as follows. We do not make any security claims
about this scheme though since that would require a separate security
analysis, and this is not in the scope of the present paper.

Private: $k\times l$ left invertible matrix $M$, with $k>l$, whose entries are
sparse polynomials from the algebra $K$.

Public: $l\times k$ right invertible matrix $L$ over $K$, such that $LM$ is
the $l\times l$ identity matrix.

Encrypting a message $m$:

  1. (1)

Convert $m$ to a vector ${\bf U}=(P_{1},\ldots,P_{l})$ of $l$ (sparse)
polynomials.

  2. (2)

Multiply the vector ${\bf U}$ by the (public) matrix $L$ on the right to get a
vector ${\bf V}=(Q_{1},\ldots,Q_{k})$ of $k$ polynomials from $K$. This vector
${\bf V}$ is the encryption of $m$.

Decryption. Multiply the vector ${\bf V}$ by the private matrix $M$ on the
right to get ${\bf V}LM={\bf U}$.

##  11\. Scrap automorphisms: scheme description

This scheme is similar in spirit to the scheme in Section 2. It is also
similar to that scheme in some key generation details, in particular in
sampling $t$-sparse polynomials. Also, converting a hash $H(m)$ to a sparse
polynomial is similar, so we chose not to duplicate Section 3.4 here.

Again, let $K=\mathbb{Z}_{q}[x_{1},\ldots,x_{n}]$ denote the algebra of
polynomials in $n$ variables over the ring $\mathbb{Z}_{q}$ of integers modulo
$q$, for an integer $q$ that is not necessarily a prime.

Let $\alpha:K\to K$ be an automorphism of the algebra $K$ given by $n$
polynomials $\alpha(x_{i})$. The signature scheme is as follows.

Public: $k$ polynomials
$y_{i_{1}}=\alpha(x_{i_{1}}),\ldots,y_{i_{k}}=\alpha(x_{i_{k}})$, where $k<n$.

Private: The inverse $\alpha^{-1}$ of the automorphism $\alpha$, given by $n$
polynomials $z_{i}=\alpha^{-1}(x_{i})$.

Signing a message $m$:

1. Apply the hash function $H$ (e.g., SHA-512) to $m$. Convert $H(m)$ to a sparse polynomial $Q$ in $y_{i_{1}},\ldots,y_{i_{k}}$ using a deterministic public procedure (see e.g. Section 3.4).

2. Apply the automorphism $\alpha^{-1}$ to the polynomial $Q$; the result is $Q(x_{i_{1}},\ldots,x_{i_{k}})$. Denote the latter polynomial by $S$.

3. The signature is the polynomial $S=\alpha^{-1}(Q)$.

Verification:

1. The verifier computes $H(m)$ and converts $H(m)$ to $Q=Q(y_{i_{1}},\ldots,y_{i_{k}})$ using a deterministic public procedure.

2. The verifier applies the public automorphism $\alpha$ to the polynomial $S=Q(x_{i_{1}},\ldots,x_{i_{k}})$.

3. The signature is accepted if and only if $\alpha(S)$ coincides with $Q=Q(y_{i_{1}},\ldots,y_{i_{k}})$.

Correctness is obvious since
$\alpha(S)=\alpha(\alpha^{-1}(Q))=\alpha(Q(x_{i_{1}},\ldots,x_{i_{k}}))=Q(y_{i_{1}},\ldots,y_{i_{k}})$.

###  11.1. Generating an automorphism $\alpha$

An automorphism $\alpha$ is generated offline, as follows. Recall that a
polynomial is $t$-sparse if it only has $t$ monomials.

(1) Choose an index $k$ at random between 1 and $n$, where $n$ is the number
of variables in our polynomial algebra.

(2) Take $x_{k}$ to $x_{k}+h(x_{1},\ldots,x_{n})$, where
$h(x_{1},\ldots,x_{n})$ is a random $t$-sparse polynomial not depending on
$x_{k}$. Fix all other variables.

(3) Apply a random permutation to the set $x_{1},\ldots,x_{n}$ of variables.

(4) Repeat steps (1) through (3) $l$ times (we suggest $l=n$).

Keep record of all “elementary” automorphisms (Steps (2) and (3)) used in the
process of building $\alpha$.

###  11.2. Computing the automorphism $\alpha^{-1}$

Having generated the automorphism $\alpha$ (see the previous subsection), we
now have to compute $\alpha^{-1}$, the private key of the signatory.

To do that, we trace back the “elementary” automorphisms (see Section 11.1)
and invert them. More specifically, if $\alpha$ is a composition
$\alpha=\varepsilon_{1}\cdots\varepsilon_{r}$ (with $\varepsilon_{1}$ applied
first), then $\alpha^{-1}=\varepsilon_{r}^{-1}\cdots\varepsilon_{1}^{-1}$.

To invert a permutation is straightforward, and if $\varepsilon_{j}$ is as in
Step (2) in the previous Section 11.1, then $\varepsilon_{j}^{-1}$ takes
$x_{k}$ to $x_{k}-h(x_{1},\ldots,x_{n})$ and fixes all other variables.

###  11.3. Sampling a $t$-sparse polynomial

The number $t$ of monomials is one of the parameters of the scheme, see
Section 11.4. Sampling a $t$-sparse polynomial can be done the same way as in
Section 11.3.

###  11.4. Suggested parameters

Again, for the integer $q$ in $\mathbb{Z}_{q}$, we suggest $q=6$.

For the number $n$ of variables, we suggest $n=32$.

For the number $k$ of published polynomials $y_{i}=\alpha(x_{i})$, we suggest
$k=16$.

For the number $t$ of monomials in $t$-sparse polynomials, we suggest $t=3$.

For the degree bound $b$ of monomials in $t$-sparse polynomials, we suggest
$b=3$.

For the number $s$ of elementary automorphisms that compose to an automorphism
$\alpha$, we suggest $s=16$.

###  11.5. Estimating the size of the key space with suggested parameters

As we have mentioned in the Introduction, the only known way to invert a
“scrap” automorphism is the “brute force” trial-and-error method. It is
therefore natural to ask for a complexity estimate of such a method. In other
words, what is the size of the key space that a brute force attacker would
have to exhaust?

First, let us estimate the number of monomials of degree $d\leq 3$ in 32
variables. By the formula for the number of combinations with repetitions,
this number is ${34\choose 3}+{33\choose 2}+{31\choose 1}=6,543$.

Thus, the number of $t$-sparse polynomials with $t=3$ is at least
${6,543\choose 3}\cdot 5^{3}>2^{42}$. (The factor $5^{3}$ appears because our
polynomial is a linear combination of 3 monomials with coefficients from
$\mathbb{Z}_{6}-\\{0\\}$.) Then the number of “elementary” automorphisms (as
in Section 11.1) is at least $2^{42}\cdot 32=2^{47}$. The number of
compositions of $s=16$ elementary automorphisms is therefore at least
$(2^{42})^{16}=2^{672}$.

Acknowledgement. We are grateful to Sasha Ushakov for useful discussions, in
particular for pointing out the forgery attack in Section 7.1.

## References

  * [1] L. Caniglia, G. Cortilnas, S. Danon, J. Heintz, T. Krick and P. Solerno, Algorithmic aspects of Suslin’s proof of Serre’s conjecture, Comput. Complexity 3 (1993), 31–55.
  * [2] A. L. Chistov, D. Yu. Grigoriev, Complexity of quantifier elimination in the theory of algebraically closed fields, in: Mathematical foundations of computer science (Prague, 1984), Lecture Notes in Comput. Sci. 176, Springer, Berlin, 1984, 17–31.
  * [3] J. Ding, T. Hodges, Cryptanalysis of an implementation scheme of the tamed transformation method cryptosystem, J. Algebra Appl. 3 (2004), 273–282.
  * [4] D. Yu. Grigoriev, M. Karpinski, M. F. Singer, The interpolation problem for k-sparse sums of eigenfunctions of operators, Adv. in Appl. Math. 12 (1991), 76–81.
  * [5] H. Lombardi, I. Yengui, Suslin’s algorithms for reduction of unimodular rows, J. Symbolic Comput. 39 (2005), 707–717.
  * [6] T. T. Moh, A public key system with signature and master key functions, Comm. Algebra 27 (1999), 2207–2222.
  * [7] NIST: Post-Quantum Cryptography: Digital Signature Schemes, https://csrc.nist.gov/csrc/media/Projects/pqc-dig-sig/documents/call-for-proposals-dig-sig-sept-2022.pdf
  * [8] D. Quillen, Projective modules over polynomial rings, Invent. Math. 36 (1976), 167–171.
  * [9] P. Shor, Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer, SIAM J. Comput. 26 (1997), 1484–1509.
  * [10] V. Shpilrain and J.-T. Yu, Polynomial automorphisms and Gröbner reductions, J. Algebra 197 (1997), 546–558.
  * [11] A. Suslin, Projective modules over polynomial rings are free, Dokl. Akad. Nauk SSSR 229 (1976), 1063–1066.