Pierre Champion

Thomas Thebaud

Gaël Le Lan

Anthony Larcher

Denis Jouvet

ON THE INVERTIBILITY OF A VOICE PRIVACY SYSTEM USING EMBEDDING ALIGNMENT

Keywords: Voice Privacy, Automatic Speaker Verification, Procrustes Analysis, Wasserstein-Procrustes

This paper explores various attack scenarios on a voice anonymization system using embeddings alignment techniques. We use Wasserstein-Procrustes (an algorithm initially designed for unsupervised translation) or Procrustes analysis to match two sets of x-vectors, before and after voice anonymization, to mimic this transformation as a rotation function. We compute the optimal rotation and compare the results of this approximation to the official Voice Privacy Challenge results. We show that a complex system like the baseline of the Voice Privacy Challenge can be approximated by a rotation, estimated using a limited set of x-vectors. This paper studies the space of solutions for voice anonymization within the specific scope of rotations. Rotations being reversible, the proposed method can recover up to 62% of the speaker identities from anonymized embeddings.

INTRODUCTION

Modern supervised deep learning algorithms require a large amount of data to be trained. To address this, service providers collect, process, and store personal data in centralized servers, raising serious concerns regarding their customer's data privacy. Recent regulations, e.g., the General Data Protection Regulation (GDPR) [START_REF]cil of 27 april 2016 on the protection of natural persons with regard to the processing of personal data and on the free movement of such data[END_REF] in the European Union, emphasize the need for service providers to ensure privacy preservation and protection of personal data. As speech data can reflect both biological and behavioral characteristics of the speaker, it is qualified as personal data [START_REF] Nautsch | The GDPR & Speech Data: Reflections of Legal and Technology Communities, First Steps Towards a Common Understanding[END_REF].

The ISO/IEC international Standard 24745 on biometric data protection [START_REF]Information technology -security techniques -biometric information protection[END_REF] defines unlinkability and non-invertibility criteria for privacy protection as follows:

• Unlinkability means that anonymized data processed in a privacy-relevant manner should not be linkable to any other set of data (anonymized or not) outside of the domain. Protected data processed in the same privacyrelevant manner must be discriminative enough to satisfy the service provider requirements, but not attackers.

• Non-invertibility means that it should be computationally infeasible 1 to obtain the clear data that led to any given anonymized data.

To achieve unlinkability, speaker anonymization [START_REF] Fang | Speaker Anonymization Using X-vector and Neural Waveform Models[END_REF][START_REF] Magariños | Reversible speaker deidentification using pre-trained transformation functions[END_REF] is performed to suppress the personally identifiable paralinguistic information from a speech utterance while maintaining the linguistic content. Recently, Fang et al. [START_REF] Fang | Speaker Anonymization Using X-vector and Neural Waveform Models[END_REF] proposed a speaker anonymization system based on the xvector paradigm and a voice conversion method. This system was used as a baseline in the first edition of the Voice privacy Challenge (VPC). The quality of anonymization is assessed using a state-of-the-art speaker verification system, which evaluates the unlinkability criteria defined in ISO/IEC 24745.

In this work, we propose to invert the VPC baseline's anonymization system using embedding alignment algorithms. In a first step, we follow the scenarios of the VPC in terms of attacker knowledge about the anonymization system [START_REF] Mohan | Evaluating voice conversion-based privacy protection against informed attackers[END_REF][START_REF] Tomashenko | Introducing the VoicePrivacy Initiative[END_REF]. The challenge assumes that the attacker has a set of clear speaker x-vectors with the corresponding anonymized x-vectors. Having this mapping allows us to approximate the anonymization function with a rotation, using a supervised Procrustes Analysis [START_REF] John | Generalized procrustes analysis[END_REF]. We propose a more restrictive scenario where the attacker does not know which clear x-vector corresponds to the anonymized x-vector. In this scenario, we use an unsupervised embedding alignment algorithm named Wasserstein-Procrustes [START_REF] Grave | Unsupervised alignment of embeddings with wasserstein procrustes[END_REF].

Once the anonymized x-vector is projected to estimate his corresponding clear x-vector, we evaluate the linkability performance between speech accessible to the attacker (enroll-ment) and speech anonymized by the service provider (trials). Invertibility is evaluated by measuring how well an attacker can invert the anonymized x-vectors of a service provider (trials). The main contributions of this paper are: (i) the approximation of a speaker anonymization system by a rotation, (ii) the use of supervised and unsupervised embedding alignment to estimate this rotation, (iii) the first (to our knowledge) invertibility attack of a speaker anonymization system.

In Section 2, we describe the VPC goals, baseline systems, data, and scenarios. Section 3 presents supervised and unsupervised alignment techniques used to perform the attack. Section 4 introduces the attack scenarios. Experimental protocol and results are detailed in Section 5 and Section 6 discusses the outcomes of this work and puts them in perspective for future research.

VOICE PRIVACY CHALLENGE

The VPC 2020 [START_REF] Tomashenko | Introducing the VoicePrivacy Initiative[END_REF] proposed an evaluation framework, dataset, and attack scenarios, which are presented in this section, to guide and facilitate the development of privacy-preserving approaches in the speech domain. The speaker anonymization system used in this work anonymizes speech segments using a x-vector-based approach [START_REF] Fang | Speaker Anonymization Using X-vector and Neural Waveform Models[END_REF]. Speaker identity and linguistic content are first extracted from an input speech utterance. Assuming that those features are disentangled, an anonymized speech waveform is generated by altering only the features that encode the speaker's identity. The anonymization system depicted in Figure 1 can be decomposed into three groups of modules. Modules from the group A extract different features from the source signal: the fundamental frequency, the phonetic features encoding articulation of speech sounds and the speaker's x-vector. The module B derives a new target identity. The x-vector from each source input speaker is compared to a pool of external x-vectors in order to select the 200 furthest vectors; 100 of them are randomly selected and averaged to create an anonymized target x-vector identity. Finally, the module C synthesizes a speech waveform from the target x-vector together with the original phonetic features and F0.

The speaker anonymization system

Speaker anonymization is achieved by selecting a private target x-vector.

Dataset

In the VPC, the evaluation dataset is built from LibriSpeech test-clean [START_REF] Panayotov | Librispeech: An asr corpus based on public domain audio books[END_REF]. Details about the number of speakers and utterances in the enrollment and trial datasets are reported in Table 1. The speech segments from 40 speakers are used to create two sets: an Enroll and a Trials set, both containing similar Female/Male ratios. Speakers from the Enroll set are all contained in the Trials sets, but their utterances are distinct between sets.

VPC attack scenarios

Speech data anonymized by the service provider is referred to as the trial dataset, and clear speech accessible to the attacker is referred to as the enrollment dataset. During the challenge, three sets of tests are performed, following the Ignorant, Lazy-Informed, and Semi-Informed attacker scenarios [START_REF] Mohan | Evaluating voice conversion-based privacy protection against informed attackers[END_REF][START_REF] Mohan | Privacy and utility of x-vector based speaker anonymization[END_REF]. In the Ignorant scenario, the attacker is unaware that speech was transformed. Thus, he performs a linkability test between the anonymized trial (denoted as Trial Anon ) and the clear, non-anonymized enrollment dataset (denoted as Enroll) using an automatic speaker verification (ASV) system with an x-vector extractor [START_REF] Snyder | X-vectors: Robust dnn embeddings for speaker recognition[END_REF] trained on the clear speech of LibriSpeech train-100 [START_REF] Tomashenko | Introducing the VoicePrivacy Initiative[END_REF][START_REF] Panayotov | Librispeech: An asr corpus based on public domain audio books[END_REF].

In the Lazy-Informed and Semi-Informed scenarios, the attacker has partial knowledge of the system and is able to anonymize the enrollment utterances using the same anonymization system but not the same target x-vector. The target x-vector chosen for each enrollment speaker differs from the target x-vector chosen for the trial speakers as the attacker does not know the randomly selected x-vectors used to generate the target identity. The difference between Lazy-Informed and Semi-Informed lies in the data used to train the x-vector extractor, for the Lazy-Informed scenario the xvector extractor is trained on clear, non-anonymized speech.

The most powerful attacker defined by the challenge is the Semi-Informed one. He has the same knowledge of the anonymization system as the Lazy-Informed, he is also able to anonymize the enrollment utterances, and additionally, he anonymizes the clear LibriSpeech train-100 dataset to retrain the x-vector extractor on anonymized data. Being computed by a retrained x-vector extractor, the x-vectors datasets used by the Semi-Informed attacker give better results.

In this paper, we provide results that follow the Lazy-Informed and Semi-Informed attacker scenarios of the VPC, but also propose additional scenarios that explore more and less constraining hypotheses for respectively unsupervised and oracle attackers.

SUPERVISED AND UNSUPERVISED ALIGNMENT ALGORITHMS

Computing the alignment of two embeddings of high dimensional real vectors is a fundamental problem in machine learning, with applications for unsupervised word and sentence translation [START_REF] Grave | Unsupervised alignment of embeddings with wasserstein procrustes[END_REF][START_REF] Biswas | Is aligning embedding spaces a challenging task? a study on heterogeneous embedding alignment methods[END_REF][START_REF] Rapp | Identifying word translations in nonparallel texts[END_REF][START_REF] Fung | Compiling bilingual lexicon entries from a non-parallel english-chinese corpus[END_REF][START_REF] Bojanowski | Unsupervised learning by predicting noise[END_REF].

Procrustes Analysis

Let A and B be two sets of N high dimensional real vectors of dimension d. We want to find the optimal rotation W ∈ R d×d that minimize the squared distance between both sets:

min W∈R d×d ||AW -B|| 2 2 (1)
For correctly matched sets A and B (the n th element of A corresponds to the n th element B, ∀n ∈ [[1, N ]]), we can directly use Procrustes analysis [START_REF] John | Generalized procrustes analysis[END_REF] to compute optimal W. This approach is well suited for supervised scenarios since it requires access to the labels of both sets. For two unlabeled sets, an unsupervised alignment algorithm is required.

Wasserstein-Procrustes

Grave et al. [START_REF] Grave | Unsupervised alignment of embeddings with wasserstein procrustes[END_REF] proposed an unsupervised algorithm to align sets of language-dependent word embeddings to perform unsupervised translation. The authors proposed a stochastic optimization, switching between minimizing the Wasserstein [START_REF] Rüschendorf | The wasserstein distance and approximation theorems[END_REF] distance between sets and finding the optimal rotation using the Procrustes analysis [START_REF] John | Generalized procrustes analysis[END_REF], to find the rotation that optimally lowers the distance between the two sets of embeddings, as well as their one-to-one mapping. In the rest of the paper, we will use this algorithm to align speaker embedding sets in unsupervised scenarios.

INVERTIBILITY ATTACK SCENARIOS

This paper explores invertibility attack scenarios (and their variations) that depend on different dataset accessibility hypotheses.

Dataset accessibility hypotheses

The datasets accessibility hypotheses are summarized in Figure 2 for the different scenarios detailed in the Section 4.2. Red boxes show data available to the attacker in a given hypothesis. Black hatched boxes show data inaccessible to the attacker in a given hypothesis. We call supervised the scenarios where labels are available and unsupervised the ones where labels are inaccessible. The scenarios where the attacker has access to the clear Trials are not realistic, but they allow to test the rotation effectiveness in the worst condition for the user. Regardless of the available data, the performances are evaluated using the Trials, Trials Anon and Enroll sets.

Scenarios

We explore invertibility attacks on the speaker anonymization system in supervised and unsupervised conditions (described in sections 4.2.1 and 4.2.2 respectively). We approximate the anonymization function in the x-vector domain by a rotation that is estimated by a supervised (Procrustes, 3.1) or unsupervised algorithm (Wasserstein-Procrustes, 3.2). Those two invertibility attacks are realistic, but we also compare them for reference to less realistic Oracle versions, where the attacker has also access to clear Trials (see Section 4.2.3).

Figure 3 (better seen in color) presents a schematic representation of the different considered attacks and or their evaluation. The Enroll sets are presented in green, the Trials sets in blue. The red arrows represent the datasets used to estimate the rotations (purple blocks). The results computed by the orange (automatic speaker verification ASV) and cyan (Top1 speaker acc.) blocks are reported in their respective lines in Table 2: the numbers next to the ASV boxes refer to the lines of Table 2 that present the corresponding results. Enroll and Trials blocks are sets of waveforms, VPC blocks refer to the x-vector-based speaker anonymization system, which takes a waveform as input, and outputs a waveform corresponding to the anonymized utterance. Purple blocks align sets of embeddings following different scenarios, which are described in this section: Procrustes (P); Wasserstein-Procrustes (WP); Procrustes oracle (PO); or Wasserstein-Procrustes oracle (WP0). Note that P, WP, PO, and WPO handle x-vectors, and that the speaker verification also relies on x-vectors. However, for better readability, the computation of the x-vectors is not explicitly shown in the figure . As an example, to obtain the results of the 12 th line of Table 2: the P matrix is computed using the Enroll and Enroll Anon sets, then applied on the Trials Anon set, the EER is computed after scoring the Enroll set against the projected P T × Trials Anon with the orange ASV block, and finally the Top 1 speaker accuracy is computed between Trials and P T × Trials Anon . For all following scenarios, once the rotation matrix W ∈ R d×d is estimated, the set of Trials Anon x-vectors is inverted using the transposed W:

Trials * W = W T × Trials Anon (2)

Supervised scenario: Procrustes

Our first scenario, the Procrustes attack, follows the rules of the VPC challenge (the hypotheses described in Section 2 are the same). It uses a rotation P, computed in a supervised manner. We apply Procrustes on the Enroll and Enroll Anon sets, knowing the one-to-one correspondence between them (thanks to the Enroll labels knowledge).

P = P rocrustes(Enroll, Enroll Anon , Enroll labels ) (3)
Then Trials Anon are inverted using equation 2 with W = P The goals of this first experiment are to measure:

• How well a rotation can approximate the VPC system in the x-vector domain by comparing the EER to the ones obtained in similar conditions for the different scenarios of the VPC.

• How many Trials Anon x-vectors can be inverted well enough to recognize their source speaker, using the reversed rotation.

Unsupervised scenario: Wasserstein-Procrustes

This second experiment explores the performance of an unsupervised algorithm for the invertibility attack. The hypothesis in this scenario presents a slight difference with the VPC ones: the attacker does not have access to the labels Enroll Labels , hence the use of the Wasserstein-Procrustes algorithm.

The optimal rotation WP can be computed using the following equation:

WP = W asserstein P rocrustes(Enroll, Enroll Anon ) (4)
The goal of this scenario is to evaluate the degradation of performance when not using Enroll Labels . Due to the VPC anonymization process, some x-vectors of the Enroll and Enroll Anon sets could be mismatched (misaligned) during the unsupervised training. Every mismatch error will contribute to degrade the alignment and lower the training and testing performances.

We also apply the variations presented in Section 4.3 to this experiment. The results are presented in lines 13 to 20 of Table 2.

Oracle scenarios

This third experiment probes the optimal performances one can get while approximating a speech anonymization system by a supervised or unsupervised rotation estimated in the xvector domain.

We suppose that an oracle has access to the Trials, Trials Anon and Trials Labels sets, meaning it can compute the best approximation possible on the evaluation data. The rotation matrix PO is computed using Procrustes:

PO = P rocrustes(Trials, Trials Anon , Trials labels ) (5)
And the rotation matrix WPO is computed using Wasserstein-Procrustes:

WPO = W asserstein P rocrustes(Trials, Trials Anon ) (6)
So the Trials Anon can be inverted with the obtained rotation matrices (same process as in equation 2).

For the oracle scenarios, the rotation is directly computed on the evaluation data. This gives the performance upper bound for the first two experiments (lines 12 and 20 of Table 2). Results are presented in lines 21 and 22 of Table 2.

The python code of the experiments is available at https:// github.com/deep-privacy/x-vector-procrustes

Experimental variations

Principal component analysis

To improve the attack performance, we extend the range of our experiments to modified x-vectors domains. We apply a dimensional reduction technique to the x-vectors sets: principal component analysis [START_REF] Jolliffe | Principal component analysis[END_REF] (PCA). Reducing the number of dimensions reduces the candidate rotations manifold, simplifying the search for the optimal one. The PCA also orders the dimensions precisely: the dimensions with the higher variance represented are placed first. This means that applying PCA on two vectors sets acts as a pre-alignment, easing the following alignment process.

For every experiment scenario using the PCA variation, we used a reduction in 70 dimensions (originally 512 for the x-vectors), and the total explained variance ratio was always above 98.0%.

Gender dependent training

As defined by the VPC evaluation rules, the unlinkability performances on Female and Male speakers are measured separately. A gender-dependent variation trains two separated rotations to improve the attack performance, one only on Female sets and the other only on Male sets.

Retrained original x-vector extractor

Corresponding to the Lazy-Informed and Semi-Informed attacker of VPC, the x-vector extractor is either trained on clear, original speech or on anonymized speech.

EXPERIMENTS

Metrics

We use two metrics to evaluate our attack on the different scenarios: Equal Error Rate (EER) and Top 1 speaker accuracy. Both metrics are computed for Female and Male speakers sets separately [START_REF] Tomashenko | Introducing the VoicePrivacy Initiative[END_REF].

For all experiments, the EER is computed by scoring the x-vectors of the reconstructed Trials Anon * set against the ones from the Enroll set, using cosine similarity. The lower the EER, the closer the reconstructed x-vectors are from the Enroll set, meaning we can find a link between the set attacked (Trials) and the set used for the attack (Enroll). A low EER would imply the capacity of the attacker to break the unlinkability aspect of the speaker anonymization system.

The Top 1 speaker accuracy is computed by comparing Trials against Trials Anon * . For each x-vector of Trials Anon * , we look for the nearest neighbor x-vector from Trials (using euclidean distance). The Top 1 speaker accuracy is the proportion of x-vectors from Trials Anon * for which the closest x-vector in Trials is from the same speaker. A high Top 1 speaker accuracy means a high success in reconstructing x-vectors close their clear counterpart and should raise concerns regarding the non-invertibility property of speaker anonymization methods.

Results

This section presents the experimental results (summarised in Table 2) for the scenarios detailed in Section 4.

We add the four scenarios presented in the voice privacy 2020 evaluation plan [START_REF] Tomashenko | Introducing the VoicePrivacy Initiative[END_REF] (see Section 2.3), for which the EER metric was recomputed using the same data but with a cosine scoring (lines 1 to 4 of Table 2). Only the Equal Error Rate is computed here because the attackers proposed in the VPC cannot inverse the speaker anonymization function.

Lines 5 to 12 explore the supervised scenario (Section 4.2.1). We can see that Procrustes gives the same attack performance as the VPC baseline attacks under the same hypothesis (similar EERs in lines 8 and 12 than 3 and 4). Regardless of the attack, the variation where the x-vector extractor is retrained on anonymized data consistently outperforms the one trained on clear data. Procrustes (line 12) achieves a Top 1 speaker accuracy of 59.8% (resp. 60.0%) for Female (resp. Male) speakers, meaning that almost six times out of 10, the anonymized speaker x-vectors can be re-identified. This raises concerns about the non-invertibility aspect of the anonymization system. The best performances are achieved by estimating a gender-dependent rotation and using PCA to reduce the x-vector dimensions (lines 8 and 12).

Lines 13 to 20 explore the unsupervised scenario (Section 4.2.2). We can see that for almost every case, Wasserstein-Procrustes gives slightly worse results than the Procrustes counterpart, as no labels are available in this scenario. We underline that the difference is usually around a few percent, so the distribution of x-vectors before and after anonymization is probably quite similar. Similar enough to get close results to when labels are available.

Lines 21 and 22 give results associated with the oracle approach (Section 4.2.3). Procrustes oracle (line 21) gives the best results among the previous experiments: 12.1% (resp. 8.7%) EER for Female (resp. Male) speakers, and 98.8% (resp. 98.0%) for the Top 1 speaker accuracy. As expected, the results are worse for Wasserstein-Procrustes oracle (line 22), with a 99.0% (resp. 98.4%) for Female (resp. Male) Top 1 speaker accuracy. Interestingly, in this scenario, with only access to both clear and anonymized x-vectors sets (no label information), the majority of x-vectors could be re-identified by the attacker.

CONCLUSION

This paper investigates various linkability and invertibility attacks on the speaker anonymization baseline of the VPC 2020. Using the challenge evaluation dataset, we approximate the anonymization function as a rotation between x-vectors domains before and after anonymization, using embedding alignment methods. Procrustes (resp. Wasserstein-Procrustes) is used to estimate the rotation in a supervised (resp. unsupervised) way. To improve performances, the attacker can compute the x-vectors thanks to an extractor retrained on anonymized utterances, estimate gender-dependent rotations, and apply PCA on the x-vectors.

Procrustes-based approaches are able to recover a large part of the mapping between clear and anonymized data; this leads to an EER which is lower than the EER calculated with the VoicePrivacy evaluation framework (line 4 and 12 of Table 2). It is also the case for the unsupervised scenario using Wasserstein-Procrustes, proving that label information is not mandatory to estimate accurate rotations. Regarding the invertibility attack, 60% Top 1 speaker accuracy is achieved in both scenarios, meaning that the inverse rotation can reidentify the majority of x-vectors. Oracle Procrustes experiment gives the upper bound for rotation approximation on the VPC baseline system: for Female and Male speakers, it could go down to 12.1% and 8.7% EER, respectively, with full access to the attacked sets.

In the unsupervised oracle scenario (e.g., full access to the unlabelled attacked sets), Wasserstein-Procrustes achieves one-to-one speaker matching between clear and anonymized counterparts with 99.0% (resp. 98.4%) accuracy for the Female (resp. Male) speakers.

The EER obtained shows that the unlinkability of a speaker anonymization system can be broken in the xvector domain by an attacker using a rotation. The Top 1 speaker accuracy leads to similar conclusions about the noninvertibility. It is of particular concern to notice that without any label information, an attacker with full access to clear and anonymized counterparts would be able to re-identify the majority of anonymized data.

Finally, these results show that there is room for improvement in current speaker anonymization systems. The unsupervised (Wasserstein-Procrustes) attack scenario seems to be an interesting approach to evaluate future anonymization methods' robustness against re-identification attacks.

Fig. 1 .

 1 Fig. 1. The Voice Privacy speaker anonymization pipeline.

Fig. 2 .

 2 Fig. 2. Schematic representation of the sets used for different scenarios. Figure best viewed in color.

Fig. 3 .

 3 Fig. 3. Schematic representation of the different attacks. Figure best viewed in color.

Table 1 .

 1 Statistics of the evaluation dataset. F and M indices refer to Female and Male speakers respectively.

	Set	Speakers Utterances Gender
	Enroll	29	438	Both
	Enroll F	16	254	Female
	Enroll M	13	184	Male
	Trials	40	1496	Both
	Trials F	20	734	Female
	Trials M	20	762	Male

Table 2 .

 2 Experimental results for the considered attack scenarios. Lines 1-4 correspond to the baseline attack scenarios of the VPC. Lines 5-12, 13-20 and 21-22 correspond to our rotation-based attack scenarios described respectively in the sections Supervised scenario (4.2.1), Unsupervised scenario (4.2.2) and Oracle scenarios (4.2.3). The variations corresponding to the columns Gender dependent, PCA and x-vector extractor are described in the section Experimental variations (4.3).

		Gender dependent	PCA	x-vector extractor	F	EER	M	Top 1 speaker Acc. F M
	1	Baseline (clear data)		original 10.3% 2.9%
	2	VPC -Ignorant		original 49.0% 42.6%
	3	VPC -Lazy-Informed		original 29.4% 29.1%
	4	VPC -Semi-Informed		retrained 17.1% 14.1%
	5			original 41.9% 30.6% 25.5%	36.6%
	6			original 40.1% 31.0% 26.6%	45.8%
	7			original 32.6% 32.7% 27.1%	40.8%
	8	Procrustes		original 25.4% 24.4% 30.5%	50.0%
	9			retrained 29.0% 23.6% 58.7%	59.2%
	10			retrained 27.0% 21.1% 54.8%	56.6%
	11			retrained 21.5% 23.1% 51.9%	57.0%
	12			retrained 14.6% 13.1% 59.8%	60.0%
	13			original 43.6% 33.4% 25.2%	22.1%
	14			original 40.7% 35.9% 26.6%	20.9%
	15			original 36.3% 35.4% 24.1%	38.6%
	16	Wasserstein		original 26.4% 25.2% 24.1%	39.4%
	17	Procrustes		retrained 31.4% 24.2% 57.2%	60.2%
	18			retrained 28.6% 24.0% 48.6%	47.1%
	19			retrained 21.6% 23.6% 48.5%	62.1%
	20			retrained 14.0% 13.2% 57.4%	61.4%
	21	Procrustes oracle		retrained 12.1% 8.7% 98.8%	98.0%
	22 Wasserstein-Procrustes oracle		retrained 13.1% 10.0% 99.0%	98.4%

Cannot be solved using an algorithm with polynomial complexity.

This work was supported in part by the French National Research Agency under project DEEP-PRIVACY (ANR-18-CE23-0018) and Région Grand Est. Experiments were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations.


Prediction of Velocity Profile in a Rectangular Shape Nozzle

Keywords: Cavitation, Multiphase, Omega, Epsilon, Turbulence ------------

Flow inside diesel injector nozzle is simulated using finite volume technique. Second order scheme is utilized to discretize the equations in order to obtain pressure in the field. It was found that K-ε turbulent model leads to more realistic results compared to K-ω turbulent model.

I. INTRODUCTION

Diesel engine is considered as a highly efficient combustion engine that plays a significant role in agriculture, transportation and industry in which numerous studies were performed in order to investigate the mentioned phenomena [START_REF] Kim | Numerical study on the cavitation phenomenon for the head drop and unsteady bubble patterns with a difference in the incidence angle of a mixed-flow pump[END_REF][START_REF] Nezamirad | Utilization of Schnerr-Sauer Cavitation Model for Simulation of Cavitation Inception and Super Cavitation[END_REF][START_REF] Sabetpour | Formation of Vapor Volume Fraction in a real size nozzle using Schnerr and Sauer approach[END_REF][START_REF] Yazdi | Effect of Pressure Difference and Needle Height on Formation of Cavitation in a real size nozzle[END_REF][START_REF] Nezamirad | Effect of Needle Height on Formation of Cavitation in a Six-Hole Diesel Injector Nozzle[END_REF][START_REF] Yazdi | Effect of Needle Height on Discharge Coefficient and Cavitation Number[END_REF][START_REF] Nezamirad | Effect of Different Diesel Fuels on Formation of the Cavitation Phenomena[END_REF][START_REF] Nezamirad | Investigation the Effect of Velocity inlet and Carrying Fluid on the Flow inside Coronary Artery[END_REF]. There have also been several studies previously that determined cavitation as an important phenomenon that should be investigated thoroughly [START_REF] Smj Zeidi | Numerical investigation of the effect of different parameters on emitted shockwave from bubble collapse in a nozzle[END_REF][START_REF] Zeidi | Investigation the effects of injection pressure and compressibility and nozzle entry in diesel injector nozzle's flow[END_REF][START_REF] Zeidi | Investigation of viscosity effect on velocity profile and cavitation formation in Diesel injector nozzle[END_REF][START_REF] Zeidi | Evaluation of the physical forces exerted on a spherical bubble inside the nozzle in a cavitating flow with an Eulerian/Lagrangian approach[END_REF][START_REF] Smj Zeidi | Effects of nozzle geometry and fuel characteristics on cavitation phenomena in injection nozzles[END_REF].

II. RESULT AND DISCUSSION

Continuity and balance of momentum can be derived as following for an incompressible flow:

𝜕𝑢 ̅ 𝑖 𝜕𝑥 𝑖 = 0 𝜕𝑢 ̅ 𝑖 𝜕𝑡 + 𝑢 ̅ 𝑗 𝜕𝑢 ̅ 𝑖 𝜕𝑥 𝑗 = - 1 𝜌 𝜕𝑃 ̅ 𝜕𝑥 𝑖 + 𝜐 𝜕 2 𝑢 ̅ 𝑖 𝜕𝑥 𝑗 𝜕𝑥 𝑗 - 𝜕 𝜕𝑥 𝑗 𝑅 𝑖𝑗
Where 𝑥 𝑖 is the position, 𝑢 ̅ 𝑖 is the mean velocity, 𝑡 is the time, 𝜐 is the kinematic viscosity, 𝑃 ̅ is the mean pressure, 𝜌 is the constant mass density, and 𝑅 𝑖𝑗 = 𝑢 ′ 𝑖 𝑢 ′ 𝑗 ̅̅̅̅̅̅̅ is the Reynolds stress tensor. Here, 𝑢 ′ 𝑖 = 𝑢 𝑖 -𝑢 𝑖 ̅ is the fluid fluctuation velocity component.

Moreover, mixture approach is used in order to calculate density in the whole domain which consist of vapor and liquid density.

ρ m =α l ρ l +α v ρ v
Where 𝑙 is refering to liquid and 𝑣 is referring to vapor.

Furthermore, viscosity of the mixture is also defined as following:

μ m =α l μ l +α v ρ v
Winklhofer [START_REF] Winklhofer | Comprehensive hydraulic and flow field documentation in model throttle experiments under cavitation conditions[END_REF] nozzle is modeled and simulated in this study. Only half of the velocity shown here is plotted. It can be inferred from figure 1 that k-ε turbulence model has a better agreement with previous experiment done by Winklhofer comparing to k-ω turbulent model. The reason for the mentioned agreement could be the differences between the law of the wall functions.

Equation ( 1) Equation ( 2)

Equation ( 3) Equation ( 4)

III. CONCLUSION

In this study, flow inside rectangular shape diesel injector nozzle is simulated using finite volume method. Continuity and momentum equations are solved using coupled scheme and the flow inside diesel injector nozzle is considered as incompressible. It was found that the velocity profile obtained from k-ε scheme has a better agreement with previous experiment done by Winklhofer comparing to k-ω turbulent model.

IV. REFERENCES

Figure 1 .

 1 Figure 1. Half velocity profile at the mid Winklhofer orifice.


Vanessa Cardin

Achim Wirth

Maziar Khosravi

Miroslav Ga

South Adriatic Recipes: estimating the vertical mixing in the deep pit

Keywords: ocean observations, vertical mixing, deep ocean, South Adriatic, time-series analysis, parameterization, climate variability

The available historical oxygen data show that the deepest part of the South Adriatic Pit remains well-ventilated despite the winter convection reaching only the upper 700 m depth. Here, we show that the evolution of the vertical temperature structure in the deep South Adriatic Pit (dSAP) below the Otranto Strait sill depth (780 m) is described well by continuous diffusion, a continuous forcing by heat fluxes at the upper boundary (Otranto Strait sill depth) and an intermittent forcing by rare (several per decade) deep convective and gravity-current events. The analysis is based on two types of data: (i) 13-year observational data time series (2006 -2019) at 750 m, 900 m, 1000 m and 1200 m depths of the temperature from the E2M3A Observatory and (ii) 55 vertical profiles (1985 -2019) in the dSAP. The analytical solution of the gravest mode of the heat equation compares well to the temperature profiles, and the numerical integration of the resulting forced heat equation compares favorably to the temporal evolution of the time-series data. The vertical mixing coefficient is obtained with three independent methods. The first is based on a best fit of the long-term evolution by the numerical diffusion-injection model to the 13-year temperature time series in the dSAP. The second is obtained by short-time (daily) turbulent fluctuations and a Prandtl mixing length approximation. The third is based on the zero and first modes of an Empirical Orthogonal Function (EOF) analysis of the time series between 2014 and 2019. All three methods are compared, and a diffusivity of approximately κ = 5 • 10 -4 m 2 s -1 is obtained. The eigenmodes of the homogeneous heat equation subject to the present boundary conditions are sine functions. It is shown that the gravest mode typically explains 99.5% of the vertical temperature variability (the first three modes typically explain 99.85%) of the vertical temperature profiles at 1 m resolution. The longest time scale of the dissipative dynamics in the dSAP, associated with the gravest mode, is found to be approximately 5 years. The first mode of the EOF analysis (85%) represents constant heating over the entire depth, and the zero mode is close to the parabolic profile predicted by the heat equation for such forcing. It is shown that the temperature structure is governed by continuous warming at the sill depth and deep convection and gravity current events play less important roles. The simple model presented here allows evaluation of the response of the temperature in the dSAP to different forcings derived from climate change scenarios, as well as feedback on the dynamics in the Adriatic and the Mediterranean Sea.

INTRODUCTION

The paradigm of the vertical temperature structure of the deep ocean is given by the "Abyssal recipes" of [START_REF] Munk | Abyssal recipes[END_REF]. It is based on a balance between an upward advective cooling and a downward diffusion of heat, which leads to an exponential profile of the temperature in the vertical direction. It explains large parts of the world's ocean temperature and salinity stratification below the thermocline. The vertical turbulent diffusion of heat needed to maintain the structure is at least two orders of magnitude larger than the molecular diffusion values of sea water, and today, it is not entirely clear where the energy that supports the vertical mixing originates from [START_REF] Munk | Abyssal recipes ii: Energetics of tidal and wind mixing. Deep-sea research. Part I[END_REF]). Tides are supposed to play an important role in the mixing process due to the breaking of internal tidal waves, generated by the interaction of the barotropic tide with the topography (see, i.e., [START_REF] Vic | Deep-ocean mixing driven by small-scale internal tides[END_REF] and de Lavergne et al. ( 2019)). Mixing in the boundary layers and generation by wind-driven horizontal turbulence are other possible sources of mixing (see [START_REF] Vallis | Atmospheric and oceanic fluid dynamics[END_REF] for a pedagogical discussion). These considerations can be applied to the Adriatic Sea (see fig. 1), whose dynamics have been widely studied in a large number of publications (see, e.g., the review by [START_REF] Cushman-Roisin | Physical oceanography of the Adriatic Sea: past, present and future[END_REF]). As far as the deep South Adriatic Pit (dSAP) is concerned, historical studies have shown that even though the winter convection very rarely goes deeper than approximately 700 m, the bottom layer is well oxygenated (see e.g., Zore-Armanda et al. (1991), [START_REF] Bensi | Thermohaline variability and mesoscale dynamics observed at the E2M3A deep-site in the South Adriatic Sea[END_REF]). Ventilation events are associated with the outbreaks of the cold North Adriatic deep water but do not occur every year. Thus, the ventilation should also be associated with some other processes, such as vertical diffusion.

More specifically, we consider the temperature structure in the area of the dSAP based on a considerable number of observations, starting from 1985 and increasing in number over the last two decades. The dSAP has a maximal depth of approximately 1230 m and is delimited at the south by the Otranto Sill (780 m deep)(see fig. 1). We focus on the temperature structure below the sill depth. It is significantly different from an exponential structure, as is shown below, and the theory put forward by [START_REF] Munk | Abyssal recipes[END_REF] has to be supplemented. The Adriatic is also interesting because the tidal motion is negligible (see, e.g., [START_REF] Cushman-Roisin | Physical oceanography of the Adriatic Sea: past, present and future[END_REF] and references therein), which simplifies the internal dynamics. An advantage of considering the dSAP is that it is easily accessible (only approximately 60 nautical miles from the coast) and has been extensively observed with considerable data coverage in space and time over the last 40 years. The present work relies on two types of data: first, in-situ ship-borne data measurements, 55 vertical temperature and salinity profiles at a 1 m resolution, collected in the dSAP from 1985 to 2019 [START_REF] Cardin | Variability of water mass properties in the last two decades in the south adriatic sea with emphasis on the period 2006-2009[END_REF]) and second, time series of high-frequency point measurements from 2006 to 2019 at 750 m, 900 m, 1000 m and 1200 m nominal depths from the mooring at the E2M3A Observatory, (http://nettuno.ogs.trieste.it/e2-m3a/). The latter are positioned in the dSAP and are particularly dedicated to the study of the long-term thermohaline and biogeochemical properties. In the present work, the diffusivity processes in the deep layer are analyzed, and thus, only measurements at 750 m (considered as the reference level for the Otranto sill depth) and deeper are taken into consideration.

It was shown by [START_REF] Querin | Saw-tooth modulation of the deep-water thermohaline properties in the southern adriatic sea[END_REF] that the dynamics in the dSAP are forced by intermittent gravity current and deep convection events that reach below the sill depth every few years and that only last from a week to a few months. The most recent events of this type occurred in March and April of 2012 and in March 2017. The former event is due to cold and very dense water from the northern Adriatic [START_REF] Bensi | Effects of winter convection on the deep layer of the southern adriatic sea in 2012[END_REF]), whereas the latter event is due to salty and warm water and has, so far, not been discussed in a shown by [START_REF] Querin | Saw-tooth modulation of the deep-water thermohaline properties in the southern adriatic sea[END_REF]. The different time scales associated with short-lasting convection/gravity current events and continuous diffusion allow separately consideration of the two regimes. In the present work, we investigate continuous mixing based on observational data using the dynamics of the heat equation.

The aim is to show that in between deep convection/gravity current events, the dynamics are governed by turbulent diffusion and that a bulk mixing coefficient can be determined to allow efficient parameterization.

The theory is presented in the next section, followed by the description of the experimental data and methodology in section 3. The results for the profiles are given in section 4.3, and those for the time series are given in section 4.4. The results are discussed in section 5, and conclusions are given in section 6.

THEORY

We consider the vertical structure of the dSAP and ignore the horizontal variations. This is justified by the dominance of horizontal (isopycnal) mixing over vertical (diapycnal) mixing [START_REF] Ledwell | Mixing of a tracer in the pycnocline[END_REF]).

Conceptually, the Munk solution is not valid in the dSAP, as there is no continuous injection in the deep waters, which justifies a continuous upward velocity. Furthermore, fig. 1 shows that the vertical gradient of the temperature vanishes towards the bottom, which indicates intense mixing without restratification effects. The data from a depth of 750 m down to the bottom have a clearly non-exponential functional form, indicating that the Munk theory of a balance between the upward advection and downward diffusion of heat does not apply. Without upward advection, only vertical turbulent diffusion remains, with no process to balance it and the only stationary solution is a constant temperature in space and time. This is not observed.

Instead, the density structure is evolving in space and time with dynamics described by the heat equation, as is demonstrated below.

Heat Equation

We are interested in the evolution of the temperature T (z, t) in the dSAP deeper than the sill depth (≈ 780 m). The underlying assumption of the present work is, that the dynamics in the dSAP are governed by the forced heat equation:

∂ t T (z, t) = κ∂ zz T (z, t) + I(z, t), (1)
where the only free parameter in the system is the vertical eddy diffusivity κ. The last term represents the intermittent intrusions in the form of deep convection or gravity-current events and is supposed to vanish over extended periods in time. The dynamics evolve between z s = -750 m, the measurement point closest to the sill depth, and the bottom, z b = -1220 m. The boundary conditions are T (z, t) = T obs (z s , t) at the sill depth and ∂ z T (z b , t) = 0 at the bottom. This shows that the model is forced continuously by a time-dependent heat flux at its upper boundary that deeply diffuses, whereas there is no heat flux at the bottom. This model is called the diffusion-intrusion model in the sequel. It is linear and therefore allows us to separately obtain analytical solutions for different types of forcing and then add them to obtain the overall solution.

The unforced heat equation, eq. ( 1), is given by I(z, t) = 0 and T obs (z s , t) = const. In the long term, the temperature relaxes to the imposed temperature at the sill depth down to the bottom, the only time independent solution. The time-dependent evolution is given by T h (z, t) = T obs (z s ) + r(t)m(z), where r(t) and m(z) are related by:

∂ t r r = κ ∂ zz m m = λ (2)
and λ does not depend on z or t, as m is not dependent on t and r is not dependent on z. (Please see Zachmanoglou and Thoe (1986) chapter 9.2 for more details on the calculations and the solution using the "separation of variables"). We further define the following: ∆T h (z, t) = T h (z, t) -T obs (z s , t). The boundary conditions of ∆T h (z, t) = 0 at the upper boundary and a vanishing temperature gradient at the lower boundary lead to the solution of the unforced heat equation of the form (z is given in meters):

∆T h (z, t) = k A k sin( π(2k -1)(z -z s ) 2(z s -z b ) ) exp(-µ(2k -1) 2 t). (3)
Note that only the wavenumbers, with odd integers 2k -1, satisfy the bottom boundary-condition and eq.

(2) imposes

µ = κ π 2(z s -z b ) 2 . ( 4
)
The characteristic e-folding decay time scale for a vertical eddy diffusivity of κ = 5 • 10 -4 m 2 s -1 (see section 4) and for the first mode, k = 1, is µ -1 = κ -1 (π/2) -2 (470m) 2 ≈ 5.7 years. For the following modes satisfying the boundary conditions, k = 2, 3 the time scales are 9 and 25 times shorter, respectively;

the modes are damped in less than a year. The higher order modes (k > 3) are damped within less than a month. The interannual dynamics are therefore dominated by the mode k = 1 (see section 4). This shows that the unforced dynamics consists of a sum of noninteracting sine modes that decay exponentially in time at a rate that is proportional to the square of their wavenumber.

V. Cardin et al.

South Adriatic Recipes

When a time-dependent forcing is applied, the dynamics of each mode are still independent. The coefficients are no longer constant in time, and we define

B k (t) = A k (t) exp(-µ(2k -1) 2 t). The first coefficient B 1 (t)
is obtained by projection on the first mode:

B 1 (t) = 4 π(z s -z b ) zs z b ∆T (z, t) sin( π(z -z s ) 2(z s -z b ) )dz (5)
The projection of the data for modes k = 2 and 3 and for data records that do not go to the total depth of 1220 m is given in appendix 1. The temporal evolution of the coefficients B k is governed by the initial conditions at t = 0, the heat flux at the sill depth and the intermittent forcing.

In the case of the intermittently forced heat equation, through deep convection or gravity currents, the temperature structure changes, and higher order modes may also have large amplitudes. This intermittent forcing is represented in eq. ( 1) through the term I(z, t). After a short (several months) intermittent event, the dynamics evolve unforced, and the amplitude of the higher modes decreases rapidly, as explained above.

Only a few months after the intermittent forcing, the dynamics are again dominated by the lower order modes.

When the continuously forced heat equation is driven by increasing temperature at a depth of 750 m at a constant rate γ = ∂ t T obs (-750m, t), the solution is:

T cf (z, t) = γt + γ 2κ (z -z b ) 2 , (6)
which can be determined by introducing the above equation in eq. ( 1) with a vanishing inhomogeneous term, I(z, t) = 0. This corresponds to the constant temperature increase of a parabolic shaped temperature profile. The temperature difference between the sill depth and the bottom of the dSAP is constant in time and is related to the temperature increase at the sill depth by:

T cf (z s , t) -T cf (z b , t) = γ 2κ (z s -z b ) 2 (7)
It is important to note that the quadratic shape is not very different from the first sine mode (see fig.

2), and when it is projected on the sine modes, the amplitudes decay quickly with the wavenumber as (2k -1) -3 . The first sine mode explains more than 92.7% of the total variance; for the first three modes, it is almost 97%. This means that from the shape of a single temperature profile at a given time t0, we cannot significantly determine whether the dynamics are freely decaying or forced by a constant-in-time vertical temperature flux. In other words, T cf (z, t 0 ) resembles T h (z, t 0 ) at a time t 0 ; however, the evolution in time for the unforced and continuously forced cases is different. In the former case, the amplitude decays exponentially, while in the latter, it stays constant.

Heat Content and Heat Flux

When a linear equation of state is supposed, the heat content, its evolution, the potential energy and its evolution as well as the increase in potential energy by mixing can be obtained independently for every mode k and are determined by the coefficients B k . If the vertical temperature profile is close to the first sine mode, the heat flux per unit area at the sill depth is

H = c p κρ∂ z T (z s ) = c p κρ π(T (z s ) -T (z b )) 2(z s -z b ) (8)
and the heat content per unit area is:

H = c p ρ 2(z s -z b ) π (T (z s ) -T (z b )). (9)
The ratio of the heat content and the heat flux is the characteristic time scale of the first mode, which is given in eq. ( 4).

EXPERIMENTAL DATA AND METHODS

In order to validate Munk's theory using experimental data we have considered the area of the South Adriatic (see fig. 1) where mixing and winter convection play an important role in homogenizing the physical and chemical properties of seawater and spreading heat along the water column. To analyze the importance of these processes, two different data sets are used: The first one consists of temperature profiles, taken during ship cruises, covering a long period of time and a large spatial area, and the second is times-series data, with high frequency point measurements along the water column by the E2M3A

Observatory. The latter is located in the dSAP and is particularly dedicated to the study of long-term thermohaline and biogeochemical properties. In this study, we analyze the diffusivity processes in the deep layer; therefore, we consider only measurements deeper than 750 m.

Profiles

The long-term variability of the thermal characteristics of the deep layer in the South Adriatic area Observatory (black square in fig. 1). Further information about these cruise operations, materials and methods can be found at http://www2.ogs.trieste.it:8585/biblioteca/.

Time Series

The data collected at the E2M3A Observatory enables the monitoring of changes that can be linked to variations in the general circulation of the Eastern Mediterranean Sea or, on a larger time scale, to the climate variability in the area, demonstrating the importance of high-frequency measurements to resolving events and rapid processes.

The Observatory has been working almost continuously since the end of 2006, and therefore, the highfrequency time series of physical and biochemical parameters measured at the site can be considered the longest almost continuous offshore data set available in the region (http://nettuno.ogs.trieste.it/e2-m3a/).

The temperature was measured using CT-CTD sensors of the SBE37 at a sampled rate of 1 hour and SBE16

with a sampling interval of 3 hours. For homogeneity, all the time series were averaged on a daily basis, and the temperature values were converted to the potential temperature. A low-pass filter PL33 with a cutoff period of 33 hours [START_REF] Flagg | 1974 mit new england shelf dynamics experiment (march, 1974)[END_REF]) was applied to remove inertial oscillations (inertial period of 17.593 hours at 42 0 N in the South Adriatic Pit) so that tidal and higher frequency fluctuations do not mask the low-frequency fluctuations in the time series driven by winds and the density field. The filter removes more than 99% of the amplitude at the semidiurnal tidal periods and more than 90% of the amplitude at the diurnal tidal periods. These subinertial non-tidal time series are considered as the starting point for the analysis. All instruments were calibrated and controlled for quality before and after each recovery and deployment, following the procedures indicated in [START_REF] Coppola | Handbook of best practices for open ocean fixed observatories[END_REF] and [START_REF] Pearlman | Evolving and sustaining ocean best practices and standards for the next decade[END_REF]. In addition, during recovery or deployment operations, a CTD cast was performed to assess the accuracy of the measured data and to identify possible instrument drift.

Empirical Orthogonal Function Analysis

To assess the temporal variability and the vertical structure of the water column in the dSAP, we applied empirical orthogonal function analysis (EOF) [START_REF] Preisendorfer | Principal component analysis in meteorology and oceanography[END_REF]) to the temperature time series. For this purpose, time series at 750 m, 900 m, 1000 m, and 1200 m were analyzed.

Unfortunately, the time series shows gaps (not simultaneous) at different levels due to instrument failures, especially in the deepest part of the mooring line in the period from 2006 to 2010 and from 2011 to 2014.

Since the EOF analysis requires data at all levels for its application, the data collected before 2014 can not used for this analysis. Thus, we considered only the time series from September 2014 to August 2019, as it was the longest period without important gaps (only short ones -no longer than 5 days -due to the Observatory service). To eliminate these gaps, the daily time series were subsequently averaged on a weekly time scale. The filtered time series allowed us to study the time scales of interest (e.g. seasonal and interannual).

RESULTS

Numerical Integration

We performed a numerical integration of the heat equation with an imposed temperature at approximately the sill depth T obs (-750m, t) (Dirichlet boundary condition) and a vanishing gradient at the bottom (1220 m) (Neumann boundary condition). For illustration, an unforced integration that starts from a linear stratification is represented in fig. 2. All values are divided by the temperature difference between 750 m and 1000 m depths to emphasize the evolution of the shape and the convergence to the first mode. Snapshots are shown at one and six months and at one, two and ten years (κ = 5 • 10 -4 m 2 s -1 ). The convergence to the first sine-mode, shown in fig. 2, with a characteristic time scale of several months, is clearly visible. In the same figure the shape of the parabolic profile for the constant forcing case is also given for comparison.

In the numerical integration, the two types of forcing, with time-dependent temperature variation at the sill depth γ and time and space dependent intermittent forcing I(z, t), can be easily combined. To reproduce the space-time dynamics of the observed data, both types of forcing are applied, and eq. ( 1) is integrated numerically.

Eddy Diffusivity

The only free parameter in the heat equation is the diffusivity. The vertical exchange of heat is at least two orders of magnitude larger than that predicted by molecular motion only, and most is performed by turbulent fluid motion. Its parameterization is often achieved by assuming a vertical eddy diffusivity. It cannot be calculated from first principles but has to be determined from observations. It can be obtained based on different data analysis methods. The first method is a best fit of the time evolution of a numerical integration described in section 4.1 to the time-series data discussed in section 3.2. More precisely, we numerically solve eq. ( 1), which is forced at the upper boundary by imposing the temperature observed at 750 m depth and supplemented by two intermittent injections in 2012 and 2017, modeling the deep convection/gravity current events that penetrate the dSAP. The value of the diffusion constant κ is then estimated so that the evolution of the model best fits the time-series data. The second method is based on Prandtl's mixing-length concept (see: Prandtl (1925), [START_REF] Taylor | Diffusion by continuous movements[END_REF][START_REF] Taylor | The present position in the theory of turbulent diffusion[END_REF] and [START_REF] Bradshaw | Possible origin of prandt's mixing-length theory[END_REF]). We first calculate a typical r.m.s. vertical scale ∆z rms = (T ) 2 /∂ z T based on the temperature fluctuations T at 900 m depth and an average temperature gradient ∂ z T at the same depth. We also determine the correlation time of the vertical motion by looking at the autocorrelation function of the temperature fluctuations

C(τ ) = T (t + τ )T (t) t / (T (t)) 2
t , which is independent of time t for statistically stationary dynamics.

We then obtain the correlation time:

τ corr = τ 0 0 C(τ )dτ
, where τ 0 is the first zero crossing of C(τ ). The vertical eddy diffusivity can then be estimated by:

κ turb = C s (∆z rms ) 2 τ corr , ( 10
)
where C s is a non-dimensional constant. The third method is based on the continuously forced dynamics, which leads to a balance between the temperature increase and the temperature difference between the sill depth and the bottom, as expressed by eq. ( 7). The three approaches for estimating the vertical eddy diffusivity are performed in sections 4.4 and4.5.

Profiles

When considering the 55 temperature profiles taken during the period from 1985 to 2019, which are presented in fig. 1, a similarity of most of them to the first sinusoidal mode depicted in fig. 2 is discernible.

For a quantitative analysis, we calculated the coefficients B k for k =1, 2 and 3, which are given in fig.

3. Only the amplitude of the first mode has a value which is significantly different from zero, varying around 0.5 o K, while the amplitudes of the second and third mode are typically more than an order of magnitude smaller and have no definite sign. We further calculated the remaining variance not explained by the first mode or the first three modes in fig. 3. For over 2/3 of the profiles, more than 99.5% of the variance can be explained by the first mode and more than 99.85% by the first three modes. Note that the 2017 event is clearly depicted by the atypical profile from 22/03/2017, given in red in fig. 1. The event also has a clear signature in the amplitude of the second and third modes (fig. 3) and in a strong increase in the variability neither explained by the first mode nor by the first three modes (fig. 3). This means that for this profile, higher order modes have an increased amplitude. In the subsequent profile only 4 months later (27/07/2017), the amplitude of the second and third mode as well as the unexplained variability has decreased to the values prior to the event. This is coherent with the fact that the decay times of modes two and three are only a few months, as explained in section 2.1. The decay of even higher modes happens at scales shorter than a month. The 2012 event, however, has no conspicuous signature in the amplitude or variability data. The profile from 31/03/2012 (given in blue in fig. 1)) shows that the dense gravity current water is just starting to arrive, and it is seen at the very bottom of the data but does not significantly impact the overall shape of the profile. No profile is available for the rest of the year, and the next profile is one year later (23/03/2013); at this time, the second, third and higher modes have faded away.

Time Series

We first numerically integrate the diffusion-intrusion model. An integration with κ fit = 5 • 10 -4 m 2 s -1 and two convection or gravity-current events in March 2012 centered at the bottom and in March 2017 centered at z = 1100 m manages to reproduce, to a reasonable degree, the observed temperature and salinity evolution, as shown in fig. 4.

Furthermore, the data at 750 m and 1200 m allows determining the amplitude of the first mode. Using this information and supposing that the higher-order modes are of low amplitude, allows to interpolate for the depths in-between. The obtained values were compared to the observed data at 900 m and 1000 m.

Again, we remark a good correspondence between the observational data and the interpolation in fig. 5. A slight deviation is observed between years 2014 and 2017.

The eddy diffusion can be estimated based on the correlation time and the vertical r.m.s displacement of a fluid parcel, as explained in section 4.2. The former is given by the integral of the correlation function from zero to the first zero crossing. The first zero crossing is at 4 days, and the integral gives τ corr = 1.6 days (not shown). The latter is calculated based on the temperature fluctuations and the mean temperature gradient (box-averaged over 30 days); its average value lies around ∆z rms = 20 m. When we impose κ turb = κ fit we find, using eq. ( 10), that C s ≈ 0.18. There is no surprise that it is smaller than one, as some of the vertical motion is due to internal waves, which do not lead to mixing if they do not break. Interestingly, the former is obtained based on turbulent fluctuations on a time scale of a few days, while the latter is obtained based on the evolution of the temperature profile over several years.

EOF analysis

The modes of the EOF analysis of the time series and the temporal evolution of EOF mode 1 are presented in fig. 6. Note that the time average is given by EOF mode 0 and that its profile resembles the parabolic profile put forward in section 2.1 for the case with a constant forcing. The dominant evolution is given by EOF mode 1 (describing 87% of the variance), which is close to a constant-in-depth temperature increase.

In this case, the temperature increase at the sill depth and the amplitude of mode 0 are related by eq. ( 7). The amplitude for mode 0 is ∆T = 0.5 o K (see fig. 6a). The warming trend is obtained by taking the product of the vertical average of mode 1 (which is close to constant see fig. 6b) and the trend of its amplitude (see fig. 7), which is close to κ turb and κ fit . We remind that the Otranto Sill depth is 780 m.

The higher order EOF modes 1, 2 and 3 closely resemble discretized versions of the analytic sine modes k=1, 2 and 3 (wavenumbers 1,3 and 5). If we suppose that they are subject to random forcing of equal strength and diffusive damping, their relative variance should decrease as 1, 1/9 and 1/25. The observed values are 1, 1/8 and 1/48, which is within discretization and statistical error (the variance of the last EOF mode is only 0.24%).

DISCUSSION

All the above results based on temperature measurements show that the action of turbulent mixing in the dSAP can be modeled well by a heat equation with a constant diffusion coefficient in space and time, a prescribed heat flux at the sill depth and a forcing term that represents rare gravity current and deep convection events. The similarity between the empirical EOF analysis results (applied only to a subset of the time series between 2014 and 2019) and a solution of the forced heat equation, with a forcing that imposes a warming trend through a heat flux at 750 m depth, confirms that the heat equation satisfactorily describes the processes in the dSAP.

When the functional form of the first mode for the temperature stratification is assumed to have an amplitude of 0.5 o K, the heat flux of the dSAP is H = 3.3Wm -2 according to eq. ( 8), and the heat content is H = 6 • 10 8 Jm -2 according to eq. ( 9) . These values characterize the heat-budget of the dSAP and can help to validate numerical models and climate dynamics.

The power per square meter necessary to perform turbulent mixing of the mass (based on temperature only) is

P = gαρκ turb ∆T. (11)
V. Cardin et al.

South Adriatic Recipes

For a ∆T = 0.5 o K, this parameter is P = 5 • 10 -4 W/m 2 . A mixing efficiency of 0.2 leads to a drain of five times more kinetic energy (see, i.e., [START_REF] Osborn | Estimates of the local rate of vertical diffusion from dissipation measurements[END_REF] and [START_REF] Peltier | Mixing efficiency in stratified shear flows[END_REF]). If the turbulent Prandtl number is unity, we obtain a characteristic vertical shear of ≈ 3 • 10 -3 s -1 .

The presented theory and observations show that from 2006 to 2019, the temperature stratification in the dSAP is determined by the continuous warming in the upper part (above the Otranto Sill depth) and the vertical turbulent diffusion. Other processes such as deep convection and gravity currents that penetrate the dSAP are only visible during a few months on two occasions. We obtain an upper (diffusive) limit for the characteristic time scale for the density structure of the dSAP, which is important as it tells us that the long-term changes in the circulation in the area are due to external forcing. The upper 5-year limit is also important for numerical modeling, as it gives a maximum spin-up time for numerical integrations. When the first mode is correctly estimated in the initial conditions of the numerical integration, a spin up of less than a year becomes sufficient, as only modes higher than k = 1 have to spin up.

The main result is that at a high accuracy, the spatio-temporal evolution of the temperature in the dSAP is described by a heat equation based on a constant (in space and time) vertical diffusivity of κ = 5 • 10 -4 m 2 s -1 . This value is similar to published values of vertical diffusivities in other regions of the world oceans. [START_REF] Munk | Abyssal recipes[END_REF], considering several stations in the interior Pacific Ocean over much larger depths, found that the dominant vertical dynamics acting on the temperature were advective diffusive, with κ = 1.3 • 10 -4 m 2 s -1 . The eddy diffusivity in deep basins of the North Aegean was estimated to lie between κ = 2 • 10 -4 to 10 • 10 -4 m 2 s -1 below a depth of 500 m ((Zervakis et al. (2003)). In the deepest part of the Calypso Deep (below 4400 m), the vertical diffusivity was found to vary around κ = (4 ± 3) • 10 -3 m 2 s -1

by [START_REF] Kontoyiannis | Hydrography, circulation, and mixing at the calypso deep (the deepest mediterranean trough) during 2006-09[END_REF].

CONCLUSIONS AND PERSPECTIVES

We showed that the theory proposed by the "Abyssal recipes" of [START_REF] Munk | Abyssal recipes[END_REF] for the world ocean is valid for the dSAP, a small oceanic basin, when forcing is added. The hypothesis of the present work is that in the dSAP disappears when forcing resides, but it represents an upper bound for the absence of deep convection events penetrating the dSAP. More precisely, when the temperature at the sill depth does not increase and no deep convection event reaches below the Otranto Sill depth within 5 years, the dSAP becomes essentially unstratified, and subsequent convection is likely to reach the bottom of the dSAP. In other words, a constant warming of the water mass above the sill depth prevents deep convection events from reaching the dSAP, and only diffusion remains; however, when there is no warming at the sill depth, a deep convection event will occur at least every five years.

The above developed theory about the dynamics of the dSAP allows a simple injection diffusion model of the dSAP to be built. This can be used to study the evolution of the dSAP under climate change scenarios, as well as to parameterize the feedback of the dSAP on the surrounding ocean dynamics. More precisely, the decomposition of modes with different characteristic time scales can be used to parameterize the response of the deep South Adriatic Pit to a specific process. For interannual variability, only the dynamics of the first mode have to be considered, while for a seasonal forcing, only the first two modes have to be taken into account, considerably reducing the complexity of the problem. The next step is to attempt to extrapolate the above findings to the more involved dynamics above the sill depth, where water masses are subject to horizontal advection and air-sea fluxes.

The performance of the theory and data analysis developed in this publication are not limited to the dSAP but can be equally applied to other enclosed pits in the world ocean, away from the direct influence of atmospheric forcing and strong currents. The temperature variation in the dSAP is significantly different from the exponential profile put forward by [START_REF] Munk | Abyssal recipes[END_REF]. As in [START_REF] Munk | Abyssal recipes ii: Energetics of tidal and wind mixing. Deep-sea research. Part I[END_REF], we have to ask the question of which process provides the power to perform the mixing in the case of the dSAP, where tides are excluded. This is one of the key open questions in physical oceanography. To answer this question, the synergy of laboratory experiments on the Coriolis platform (Grenoble) of idealized numerical simulations and of fine resolution (in time-and-space) velocity measurements in the dSAP is currently being employed.

The proposed empirical model for the vertical, one-dimensional, structure of the temperature profile, its evolution in time, as well as the determination of the vertical mixing coefficients and the forcing, relies on the synergy between mooring data with dense temporal coverage and a substantial number of profiles, providing fine resolution in the vertical direction. This dataset has been collected during the last four decades, at a single oceanic site and shows the importance of a continuous effort of collecting data over several generations of funding programs.

PROJECTION

The discretized version of the projection to obtain the amplitude of the modes is given for the first three modes, it is defined iteratively by: If the record is over the total length λ k = 1. If the record stops at length z B < 1219m, so does the summation in eqs. ( 12) -( 14) and:

B 1 (t)
λ k = z B -z s 469m - sin(2(k -1)π z B -zs 469m ) 2(k -1)π • 469m -1 (15)

Figure 1 .

 1 Figure 1. Left: Region of the South Adriatic Pit (SAP). The red rectangle shows the area where all the profiles originated. After 2010, data were collected within the black rectangle. The yellow dot indicates the position of the E2M3A observatory of the thermohaline time series. Right: All 55 profiles of the temperature dependence with the depth. The two profiles during the deep convection/gravity current events are emphasized in blue and red.

  is determined based on temperature data collected during several oceanographic campaigns between 1985 and 2019. The cruises were undertaken as part of Italian and European projects, and other datasets were kindly provided by the NATO UNDERSEA RESEARCH CENTER, PANGEA Data Repository (https://www.pangaea.de/) and MEDATALAS Dataset (https://nodc.inogs.it/nodc/). All the profiles were processed following the methodology described by[START_REF] Cardin | Variability of water mass properties in the last two decades in the south adriatic sea with emphasis on the period 2006-2009[END_REF] and applied for the data between 1990 and 2010, i. e. the oceanographic condition during a specific campaign is represented by a profile obtained by averaging all the CTD casts deeper than 1000 m measured in the area between 41.25 o and 42.25 o N and 17.25 o and 18.25 o E in the SAP (fig.1). For the studied period, we consider data from 55 temperature profiles deeper than 750 m. The potential temperature (hereafter referred to as temperature) is calculated from the original profiles with the surface as the reference depth. Profiles after 2010 were collected in a narrower area of 0.5 o x0.5 o degrees Lat/Lon located in the dSAP and near the E2M3A

For

  this study, the time series considered along the mooring line come from 4 instruments positioned at 750 m, 900 m, 1000 m and 1200 m(Cardin and Bensi (2014),Cardin et al. (2014),[START_REF] Cardin | E2m3a-2013-2015-time-seriessouthadriatic[END_REF] and[START_REF] Cardin | E2m3a-2015-2017-time-seriessouthadriatic[END_REF]). To simplify the analysis, the depths of the series are nominal. The available data used are from November 2006 to October 2019, with a gap between September 2010 and May 2011. The time series at 900 m started in November 2013. An exhaustive description of the site, data availability and data analysis can be found in[START_REF] Bensi | Thermohaline variability and mesoscale dynamics observed at the E2M3A deep-site in the South Adriatic Sea[END_REF],[START_REF] Bensi | Thermohaline variability and mesoscale dynamics observed at the deep-ocean observatory e2m3a in the southern adriatic sea. The Mediterranean Sea: Temporal Variability and Spatial Patterns[END_REF] and at the Observatory webpage (http://nettuno.ogs.trieste.it/e2-m3a/).

Figure 2 .

 2 Figure 2. Temperature profiles from the numerical model given by the unforced heat equation (1) (all values are divided by the temperature difference between 750 m and 1000 m depths) for the initial condition (black line), after one month (red), six months (green), one year (blue), two years (yellow) and ten years (brown) of the integration time (κ = 5 • 10 -4 m 2 s -1 ). The black dashed line is the analytic solution of the first mode, to which the dynamics converge. The red dashed line is the parabolic profile of the continuous forcing.

Figure 3 .

 3 Figure 3. Upper row: evolution of the amplitudes of the sine modes B 1 (black), B 2 (red) and B 3 (green). The vertical scale gives the temperature difference between the 750 m time series and the maximum profile depth (1219 m) of the corresponding mode. Lower row: part of variability not explained by B 1 (black), B 1 + B 2 (red) or B 1 + B 2 + B 3 (green) (left) and zoom (right). Figures on the right show the same data magnified.

Figure 4 .

 4 Figure 4. The observational data (thin lines) and the corresponding solutions of the injection-diffusion model (thick lines) for different depths (black 750 m, red 900 m, green 1000 m, blue 1200 m) for temperature.

Figure 5 .

 5 Figure 5. Observed data at 900 m (red) and 1000 m (light green) and their interpolations 900 m (orange) and 1000 m (dark green). The right figure is a magnification of the data starting in 2017.

Figure 6 .

 6 Figure 6. Temporal average of data, EOF mode 0 (a) and the EOF modes 1 to 4 with their variances given in the legend (b). Temporal evolution of the amplitude of EOF mode 1 (c).

  the dynamics in the dSAP are governed by turbulent diffusion and forced by a continuous heat flux at the upper boundary and sporadic deep convection/gravity current events. A hierarchy of models, all based on the heat equation, is used. The simplest are based on the gravest or the three gravest modes of the analytically solved unforced heat equation. It is shown that the heat equation forced by a constant heat flux at the upper boundary has a functional form close to the gravest mode. The most evolved model includes a time-dependent forcing at the upper boundary as well as intermittent convection events that are solved numerically.The models successfully describe the profiles their spatio-temporal evolution. The functional form of the profile is described by the gravest mode of the analytical solution of the homogeneous heat equation, while the time series can be fitted well by the numerical solution of the forced heat equation. Two major deep convection/gravity current events are identified in the time series in 2012 and 2017 between 2007 and 2019 (during the winter/spring season). The longest dissipative time scale in the dSAP, associated with the gravest mode, is approximately 5 years. This does not only indicate that within this period, stratification

This is a provisional file, not the final typeset article

Frontiers

ACKNOWLEDGMENTS

Part of this work was performed when AW visited OGS and VC visited LEGI. We thank the OGS ExO and TECDEV Group so as the Captain and crew of the RV OGS Explora and RV Laura Bassi for the valuable support and work during the maintenance of the E2M3A Observatory.

FUNDING

This work was funded by Labex OASUG@2020 (Investissement d'avenir -ANR10 LABX56) and OGS internal funding.

DATA AVAILABILITY STATEMENT

Data sets are referenced in the references.

CONFLICT OF INTEREST STATEMENT

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

AUTHOR CONTRIBUTIONS

All authors have contributed to the paper. The numerical projection and integration was performed by AW, the EOF analysis by MK.

Zachmanoglou, E. C. and Thoe, D. W. (1986). Introduction to partial differential equations with applications (Courier Corporation) Zervakis, V., Krasakopoulou, E., Georgopoulos, D., and Souvermezoglou, E. (2003). Vertical diffusion and oxygen consumption during stagnation periods in the deep north aegean. Deep Sea Research Part I: Oceanographic Research Papers 50, 53-71 Zore-Armanda, M., Bone, M., Dadic, V., Morovic, M., Ratkovic, D., Stojanoski, L., et al. (1991).

Hydrographic properties of the adriatic sea in the period from 1971 through 1983. Acta adriatica. 1991.