Christina Boura
email: christina.boura@uvsq.fr

Anne Canteaut
email: anne.canteaut@inria.fr

Daniel Coggia
email: daniel.coggia@inria.fr

A General Proof Framework for Recent AES Distinguishers

Keywords: AES, Distinguisher, Subspace Trail Cryptanalysis

In this paper, a new framework is developed for proving and adapting the recently proposed multiple-of-8 property and mixture-differential distinguishers. The above properties are formulated as immediate consequences of an equivalence relation on the input pairs, under which the difference at the output of the round function is invariant. This approach provides a further understanding of these newly developed distinguishers. For example, it clearly shows that the branch number of the linear layer does not influence the validity of the property, on the contrary of what was previously believed. We further provide an extension of the mixture-differential distinguishers and multiple-of-8 property to any SPN and to a larger class of subspaces. These adapted properties can then be exhibited in a systematic way for other ciphers than the AES. We illustrate this with the examples of Midori, Klein, LED and Skinny.

Introduction

The Advanced Encryption Standard (AES) is a block cipher designed by Daemen and Rijmen in 1997 and standardised by the NIST in 2001 [START_REF]Advanced Encryption Standard (AES)[END_REF]. It is since then the most used and the most analysed symmetric primitive worldwide. Since its submission to the NIST competition, many different cryptanalytic techniques, including among others, integral attacks [FKL + 01], or sophisticated meet-in-the-middle attacks [START_REF] Demirci | A meet-in-the-middle attack on 8-round AES[END_REF][START_REF] Derbez | Exhausting Demirci-Selçuk meetin-the-middle attacks against reduced-round AES[END_REF][START_REF] Derbez | Automatic search of meet-in-themiddle and impossible differential attacks[END_REF] have been developed and applied to the AES, exploiting different properties of the algorithm. Up to now, all the developed techniques, at least in the single-key model, have only managed to break reduced-round versions of the standard1 .

Most of the attacks against block ciphers are based on the existence of a distinguisher, that is a non-random property that permits to distinguish within reasonable time and by using reasonable data and memory resources, a reduced-round version of the cipher instantiated with a random secret key from a random permutation. Until recently, all known distinguishers of the AES in the single-key model could reach at most 4 rounds. However, since 2016, the first 5-round AES-distinguishers appeared [SLG + 16, GRR17, RBH17, Gra18] and this topic has become again a subject of broad and current interest. The importance of these distinguishers is that they exhibit new, unexplored properties of the AES. They led to improved attacks on reduced-round versions of the cipher, like the attack on 5 rounds described in [BDK + 18] based on the distinguisher exhibited in [START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF].

The main breakthrough in these attacks is the identification by Grassi, Rechberger and Rønjom [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF] of the following property of the AES round function R. There exist two well-chosen linear subspaces V and W of F 128 2 satisfying the following property: for any coset of V , (c + V ), the number of distinct pairs of elements x, x , x = x in (c + V ) such that R(x) and R(x ) belong to the same coset of W is always divisible by 8. This behaviour, known as the multiple-of-8 property, is then combined with two 2-round deterministic subspace trails to form a 5-round distinguisher. However, despite the theoretical interest of this distinguisher, because of the nature of the particular subspace V used, it could not be exploited directly for mounting a key-recovery attack. For this reason, Grassi presented in [START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF] new 4-round distinguishers, that exploit a property appearing in the proof of the multiple-of-8 property but that is expressed in a more convenient way and facilitates key-recovery attacks. These new distinguishers were given the name of mixture-differential distinguishers.

The proof of the multiple-of-8 property given in [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF] is divided into many special cases and each case needs to be proved separately. An adaptation of the initial property, considering slightly more general input subspaces is also provided in the same paper and it requires a similar case-by-case proof. The disadvantage of these redundant proofs is that it is not clear from them what are the characteristics and the properties of the inner components of the AES that have an influence on the multiple-of-8 behaviour. Furthermore, the question whether this kind of property is proper to the AES or whether it can be adapted to other ciphers is unclear and the original proofs do not provide hints to answer this question. The same questions can be asked for the mixture-differential distinguishers.

Then, the aim of our paper is to provide a general formulation of the mixture-differential distinguisher and of the multiple-of-8 property which can be applied in a systematic way to any cipher following the SPN construction. It then avoids all these redundant proofs which were previously required for each new occurrence of these properties. Also, our result precisely identifies the conditions to be satisfied for the property to hold. Most notably, it shows that these distinguishers apply to a more general class of SPN than the ones mentioned in [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF][START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF].

Our contributions. We show that the mixture-differential distinguishers, and by extension the multiple-of-8 property, revealed in [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF][START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF] are direct consequences of the fact that the difference between two outputs of the AES round function, R(p 0 ) + R(p 1 ), is invariant under an equivalence relation between the plaintext pairs. The definition of this equivalence relation leads to a simple and compact proof of the distinguishing properties exhibited in [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF][START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF]. Also it clarifies which parts of the AES have an influence on the property. Most notably, the validity of the property does not depend on the branch number of MixColumns, on the contrary of what was previously believed in [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF][START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF]. This disproves for instance the statement in the abstract of [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF], which points out that "this new structural property [...] is independent of the details of the MixColumns matrix (with the exception that the branch number must be maximal)". We show here that the property holds even if MixColumns has a lower branch number. Also, we describe the form of the linear subspaces that can be used in the distinguisher in place of the mixed spaces M I in the original article. This permits to adapt the property and makes it directly applicable to other ciphers. As an illustration, we show that the same kind of property holds for the block ciphers Midori [BBI + 15], Klein [START_REF] Gong | KLEIN: A new family of lightweight block ciphers[END_REF], LED [START_REF] Guo | The LED block cipher[END_REF] and

Skinny [BJK + 16].
Organisation of the paper. The rest of the paper is organised as follows. A description of the AES and basic definitions on subspace trails are recalled in Section 2. Section 3 describes the multiple-of-8 property and the mixture-differential distinguishers presented in [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF][START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF]. Section 4 then exhibits an equivalence relation between pairs of AES states, which leads to a new proof of the above properties. It also discusses the influence of the branch number of MixColumns on these results. Then, an adaptation to any SPN and other linear subspaces is presented in Section 5. Finally, applications on different ciphers are provided in Section 6.

Preliminaries

We start by providing a brief description of the AES and introduce in parallel the notation that we will use. We work throughout this paper with finite fields of characteristic 2 that are fields containing 2 d elements, seen as extensions of F 2 .

Description of the AES

The Advanced Encryption Standard [START_REF]Advanced Encryption Standard (AES)[END_REF] is a Substitution-Permutation network operating on 128-bit plaintexts. The master-key size can be 128, 192 or 256 bits and the round-key size is 128 bits. The number of rounds N r is respectively 10, 12 or 14. The internal state is represented as a 4 × 4 matrix over the field F 2 8 , called the state array. An AES round R is the composition K • L • S where:

• S is the SubBytes operation applying the same invertible S-box to each F 2 8 -entry of the state array.

• L = MC • SR is the linear layer. SR, is the ShiftRows operation consisting of a cyclic shift of each row to the left and MC, which stands for MixColumns, consists in the left multiplication of the state array by a 4 × 4 constant matrix over F 2 8 . This matrix, denoted M MC , is defined by

M MC =     2 3 1 1 1 2 3 1 1 1 2 3 3 1 1 2     ,
where F 2 8 is identified with F 8 2 based on the irreducible polynomial X 8 + X 4 + X 3 + X + 1, and the elements of F 8 2 are represented as integers in {0, . . . , 2 8 -1}.

• K is the AddRoundKey operation bitwise adding to the state array a 128-bit round-key derived from the master key.

An important notion for measuring the quality of the linear layer of a block cipher is the so-called differential branch number. This notion, that will play a role later in our analysis is introduced below.

Definition 1 (Differential branch number [Dae95]). Let M : (F d 2 ) N → (F d 2 ) N be an F 2 -linear function. The differential branch number of M with respect to F d 2 is min x∈(F d 2 ) N \{0} (wt d (x) + wt d (M x))
where

wt d (x) = #{i ∈ {0, . . . , N -1} | x i = 0} when x = (x 0 , . . . , x N -1 ) ∈ (F d 2 ) N .
Notably, AES MixColumns MC has branch number 5 with respect to F 8 2 .

Subspace trails for AES

We now recall the basic notation on subspace-trail cryptanalysis of AES, first introduced in [START_REF] Grassi | Subspace trail cryptanalysis and its applications to AES[END_REF].

Let d be the degree of the extension over F 2 on which the S-box operates and K = F 2 d . For the AES, d = 8, but we rather present all results for an arbitrary d since simulations are often performed on small scale variants of the AES with a 4-bit S-box (see e.g. Section 5.2 in [START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF]). Let M 4 (K) denote the set of all 4 × 4-matrices over K and let (e i,j ) i,j∈{0,...,3} be the canonical basis of M 4 (K):

j ↓ e i,j =    0 • • • 0 . . . 1 . . . 0 • • • 0    ← i
In the following, vect K (v 0 , . . . , v k-1 ) denotes the linear space formed by all linear combinations with coefficients in K of the vectors v 0 , . . . , v k-1 ∈ M 4 (K). As in [START_REF] Grassi | Subspace trail cryptanalysis and its applications to AES[END_REF][START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF], we define the following subspaces of M 4 (K) for i ∈ {0, . . . , 3}, with indices computed modulo 4.

The column spaces

: C i = vect K (e 0,i , e 1,i , e 2,i , e 3,i ), The diagonal spaces :

D i = vect K (e 0,i , e 1,i+1 , e 2,i+2 , e 3,i+3 ) = SR -1 (C i ), The anti-diagonal spaces : ID i = vect K (e 0,i , e 1,i-1 , e 2,i-2 , e 3,i-3 ) = SR(C i ), The mixed spaces : M i = MC(ID i ). For example, if x 0 , x 1 , x 2 , x 3 ∈ K,    
x 0 0 0 0 x 1 0 0 0 x 2 0 0 0

x 3 0 0 0     ∈ C 0 ,     x 0 0 0 0 0 x 1 0 0 0 0 x 2 0 0 0 0 x 3     ∈ D 0 ,     x 0 0 0 0 0 0 0 x 1 0 0 x 2 0 0 x 3 0 0     ∈ ID 0 ,     2 • x 0 x 1 x 2 3 • x 3 x 0 x 1 3 • x 2 2 • x 3 x 0 3 • x 1 2 • x 2 x 3 3 • x 0 2 • x 1 x 2 x 3     ∈ M 0 .
If I ⊆ {0, 1, 2, 3}, we then define :

C I = i∈I C i , D I = i∈I D i , ID I = i∈I ID i , M I = i∈I M i .
Now that we have linear subspaces of M 4 (K), we define their cosets as affine subspaces of M 4 (K). More precisely, a coset of the linear subspace V ⊆ M 4 (K) is a set of the form

V + a = {v + a | v ∈ V } where a ∈ M 4 (K)
. Moreover, we are going to pay attention to subspaces that satisfy a specific property defined below.

Definition 2 (Subspace trail [START_REF] Grassi | Subspace trail cryptanalysis and its applications to AES[END_REF]). Let F : K N → K N be any map. Two linear subspaces U, V ⊆ K N form a (one-round) F-subspace trail if

∀a ∈ K N , ∃b ∈ K N : F(U + a) ⊆ V + b, (1)
which is denoted by

U F ⇒ V . The negation is denoted by U F ⇒ V . An (r + 1)-tuple of subspaces (U 0 , . . . , U r ) is called a subspace trail (over r rounds) if ∀i ∈ {0, . . . , r -1}, U i F ⇒ U i+1 .
For example, we have the trivial subspace trails {0} F ⇒ {0} and U F ⇒ K N . In this paper, we only consider exact subspace trails, i.e., for which equality holds in (1).

Distinguishers based on subspace trails

We describe in this section three distinguishers based on subspace trail cryptanalysis in order to have a complete understanding of the context in which our contribution lies. The first one, describes the multiple-of-8 property presented by Grassi, Rechberger and Rønjom at Eurocrypt 2017 [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF], while the other two are based on the mixture differential property and were published recently by Grassi in [START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF].

The distinguisher from [GRR17]

We first describe the distinguisher presented in [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF]. We begin with this easy to verify lemma describing a two-round subspace trail for the AES.

Lemma 1 ([GRR16]). Let I ⊆ {0, 1, 2, 3}, then D I R ⇒ C I R ⇒ M I .
Now comes a more subtle lemma which is the keystone of Theorem 1. In the whole paper, as in [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF], we always consider unordered pairs of elements and denote them as pair sets, i.e. {a, b}.

Lemma 2 ([GRR17]). Let a ∈ M 4 (K), i ∈ {0, 1, 2, 3}, J ⊆ {0, 1, 2, 3}. We define n = #{ {p 0 , p 1 } with p 0 , p 1 ∈ M i + a | R(p 0 ) + R(p 1 ) ∈ D J }.
Then n ≡ 0 mod 8. A proof of Lemma 2 is given in Section 6 of [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF] but we provide in Section 4 a much more compact proof of this same result. A direct consequence of Lemmas 1 and 2 is the following theorem.

Theorem 1 ([GRR17]). Let a ∈ M 4 (K), i ∈ {0, 1, 2, 3}, J ⊆ {0, 1, 2, 3}. We define n = #{ {p 0 , p 1 } with p 0 , p 1 ∈ D i + a | R 5 (p 0 ) + R 5 (p 1 ) ∈ M J }.
Then n ≡ 0 mod 8.

Theorem 1 directly provides a distinguisher for five rounds of the AES independent of the secret key. Indeed, given an oracle simulating either five rounds of the AES, either a random permutation, one can compute the number n from Theorem 1 with only the 2 32 plaintexts belonging to the same coset of D i . This distinguisher is fully described in [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF].

Exploiting the above distinguisher for mounting a key-recovery attack on more rounds revealed to be however a difficult task because of the form of the output subspace M J . Indeed, as M J affects the whole AES-state, a key-recovery attack requires the guess of the entire subkey in the last round. For this reason, Grassi presented in [START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF] new distinguishers that exploit similar properties but have a description that is more adapted to a key-recovery attack. The counterpart is that these distinguishers cover one round fewer than in [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF].

The two distinguishers from [Gra18]

In the following, we use the notation from [START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF]: for a given basis (g 0 , . . . , g k-1 ) and a given element a ∈ M 4 (K), a vector p in the affine subspace a + vect K (g 0 , . . . , g k-1 ) defined by

p = a + k-1 i=0 x i g i with x i ∈ K is denoted by p ≡ (x 0 , . . . , x k-1 ). Theorem 2 ([Gra18]). Let a ∈ M 4 (K), i ∈ {0, 1, 2, 3}, J ⊆ {0, 1, 2, 3}
and let p 0 , p 1 , q 0 , q 1 ∈ vect K (e 0,i , e 1,i ) + a. Suppose that 1. p 0 ≡ (x 0 , x 1 ), p 1 ≡ (y 0 , y 1 ), 2. q 0 ≡ (x 0 , y 1 ), q 1 ≡ (y 0 , x 1 ).

Then R 4 (p 0 ) + R 4 (p 1 ) ∈ M J if and only if R 4 (q 0 ) + R 4 (q 1 ) ∈ M J .
We will provide an alternative proof of Theorem 2 in Section 5. The following variant of the distinguisher, involving another input subspace, is also exhibited in [START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF].

Theorem 3 ([Gra18]). Let a ∈ M 4 (K), i ∈ {0, 1, 2, 3}, J ⊆ {0, 1, 2, 3} and let p 0 , p 1 ∈ C i + a. Suppose that p 0 ≡ (x 0 , x 1 , x 2 , x 3 ) and that p 1 ≡ (y 0 , y 1 , y 2 , y 3 ). Then, R 4 (p 0 ) + R 4 (p 1 ) ∈ M J if and only if R 4 (q 0 ) + R 4 (q 1 ) ∈ M J ,
for all sets of plaintexts {q 0 , q 1 } with q 0 , q 1 ∈ C i + a of the following form:

1. q 0 ≡ (y 0 , x 1 , x 2 , x 3 ), q 1 ≡ (x 0 , y 1 , y 2 , y 3 ) 2. q 0 ≡ (x 0 , y 1 , x 2 , x 3 ), q 1 ≡ (y 0 , x 1 , y 2 , y 3 )
3. q 0 ≡ (x 0 , x 1 , y 2 , x 3 ), q 1 ≡ (y 0 , y 1 , x 2 , y 3 ) 4. q 0 ≡ (x 0 , x 1 , x 2 , y 3 ), q 1 ≡ (y 0 , y 1 , y 2 , x 3 ) 5. q 0 ≡ (x 0 , x 1 , y 2 , y 3 ), q 1 ≡ (y 0 , y 1 , x 2 , x 3 ) 6. q 0 ≡ (x 0 , y 1 , x 2 , y 3 ), q 1 ≡ (y 0 , x 1 , y 2 , x 3 ) 7. q 0 ≡ (x 0 , y 1 , y 2 , x 3 ), q 1 ≡ (y 0 , x 1 , x 2 , y 3 )

We will provide a proof of Theorem 3 in Section 4.

A more concise and general proof

This section is dedicated to our proof of Lemma 2 from [GRR17] (Lemma 2 here). This new proof is a much more concise version of the case-by-case proof given in the original paper. To be more precise, instead of proving Lemma 2, we prove directly a more general variant. This generalisation is present in the original paper but its proof is only sketched in Appendix A of [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF]. Indeed, the proof framework of Lemma 2 in [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF] does not allow a compact proof of this generalisation.

Our approach for proving Lemma 2 can be divided into three steps:

• there exists an equivalence relation between pairs of elements in a certain subspace of M 4 (K) (Definition 4);

• some function on those pairs derived from the round function is invariant under this equivalence relation (Theorem 4);

• the cardinality of the equivalence classes is always a multiple of 8 (Proposition 1).

In the following, we fix a ∈ M 4 (K), I ⊆ {0, 1, 2, 3} and J ⊆ {0, 1, 2, 3}. Here it might help to remind that

ID I = vect K (e k,i-k | k ∈ {0, 1, 2, 3}, i ∈ I) = vect K (e i-k,k | k ∈ {0, 1, 2, 3}, i ∈ I), M I = MC(ID I ) = vect K (MC(e i-k,k ) | k ∈ {0, 1, 2, 3}, i ∈ I).

An equivalence relation between pairs of states

Definition 3 (Information set). Let {p 0 , p 1 } be a set of elements in M I + a, written as

p 0 = 3 k=0 i∈I p 0 i,k MC(e i-k,k ) + a and p 1 = 3 k=0 i∈I p 1 i,k MC(e i-k,k ) + a
for some (uniquely defined)

p 0 i,k , p 1 i,k ∈ K, i ∈ I, 0 ≤ k ≤ 3. The information set K of {p 0 , p 1 } is defined as K = {k ∈ {0, 1, 2, 3} | ∃i ∈ I : p 0 i,k = p 1 i,k } .
Definition 4 (Equivalence relation). Let P = {p 0 , p 1 } and Q = {q 0 , q 1 } with p 0 , p 1 , q 0 , q 1 ∈ M I + a. We say that P ∼ Q if:

• {p 0 , p 1 } and {q 0 , q 1 } have the same information set K.

• ∀k ∈ K, ∃b ∈ {0, 1} : ∀i ∈ I, q 0 i,k = p b i,k and q 1 i,k = p 1-b i,k .
Clearly, ∼ is an equivalence relation on unordered pairs of M I + a.

Example 1. The following two sets {p 0 , p 1 } and {q 0 , q 1 }, with p 0 , p 1 , q 0 , q 1 ∈ M 0 are equivalent.

{p 0 , p 1 } =            2 • x 0 x 1 z 2 3 • z 3 x 0 x 1 3 • z 2 2 • z 3 x 0 3 • x 1 2 • z 2 z 3 3 • x 0 2 • x 1 z 2 z 3     ,     2 • y 0 y 1 z 2 3 • z 3 y 0 y 1 3 • z 2 2 • z 3 y 0 3 • y 1 2 • z 2 z 3 3 • y 0 2 • y 1 z 2 z 3            ∼ {q 0 , q 1 } =            2 • x 0 y 1 w 2 3 • w 3 x 0 y 1 3 • w 2 2 • w 3 x 0 3 • y 1 2 • w 2 w 3 3 • x 0 2 • y 1 w 2 w 3     ,     2 • y 0 x 1 w 2 3 • w 3 y 0 x 1 3 • w 2 2 • w 3 y 0 3 • x 1 2 • w 2 w 3 3 • y 0 2 • x 1 w 2 w 3           
The information set of these pairs has cardinality |K| = 2. Now that we have the right definitions, we can state and prove the following key theorem that has also other applications than just proving Lemma 2. It is worth noticing that, in the original proof of Lemma 2, the authors split the proof procedure in several cases, each case corresponding to a different size of the information set. In our approach we assemble all cases together by using the above introduced equivalence relation.

Theorem 4. For any a ∈ M 4 (K), the function ∆ operating on unordered pairs of elements in M I + a and defined by

∆ : {p 0 , p 1 } -→ R(p 0 ) + R(p 1 )
is constant over the equivalence classes for ∼.

Proof. Let P = {p 0 , p 1 } and Q = {q 0 , q 1 } with p 0 , p 1 , q 0 , q 1 ∈ M I + a such that P ∼ Q. We write as in Definition 3

p 0 = 3 k=0 i∈I p 0 i,k MC(e i-k,k ) + a and p 1 = 3 k=0 i∈I p 1 i,k MC(e i-k,k ) + a.
We also write the MixColumns matrix M MC = (m ,k ) 0≤ ,k≤3 . Hence

p 0 = k, i∈I p 0 i,k m ,i-k e ,k + a = k, i∈I p 0 i,k m ,i-k + a ,k e ,k .
Then

S(p 0 ) = k, S-box i∈I p 0 i,k m ,i-k + a ,k e ,k and
S(p 0 ) + S(p 1 ) = k, S-box i∈I p 0 i,k m ,i-k + a ,k + S-box i∈I p 1 i,k m ,i-k + a ,k e ,k .
(2)

It is now clear with Definition 4 and Equation (2) that S(p 0 ) + S(p 1 ) and S(q 0 ) + S(q 1 ) are equal in M 4 (K). Indeed, with K the information set of P and Q,

S(q 0 ) + S(q 1 ) = k, ∈{0,1,2,3} S-box i∈I q 0 i,k m ,i-k + a ,k + S-box i∈I q 1 i,k m ,i-k + a ,k e ,k = ∈{0,1,2,3}, k∈K S-box i∈I q 0 i,k m ,i-k + a ,k + S-box i∈I q 1 i,k m ,i-k + a ,k e ,k = ∈{0,1,2,3}, k∈K S-box i∈I q b(k) i,k m ,i-k + a ,k + S-box i∈I q 1-b(k) i,k m ,i-k + a ,k e ,k = ∈{0,1,2,3}, k∈K S-box i∈I p 0 i,k m ,i-k + a ,k + S-box i∈I p 1 i,k m ,i-k + a ,k e ,k = S(p 0 ) + S(p 1 ) .
Therefore,

∆(P ) = R(p 0 ) + R(p 1 ) = K • L • S(p 0 ) + K • L • S(p 1 )
= L(S(p 0 ) + S(p 1 )) by characteristic 2 and linearity of L;

= L(S(q 0 ) + S(q 1 )) by our previous observation; = ∆(Q).

We would like to adapt in Section 5 the previous theorem to any SPN structure and any linear space. For this, it is important to identify the key argument that makes the proof work. Looking at the proof carefully, it appears that it relies on the fact that the coordinates p i,k of the elements in M I can be decomposed into disjoint sets, in such a way that each individual S-box involves only one coordinate subset. Here, the four subsets correspond to the coordinates p i,k sharing the same index k. Indeed, the S-box at Row and Column k involves all p i,k , i ∈ I. This particular structure then makes possible to exchange between the two pairs {p 0 , p 1 } and {q 0 , q 1 } the coordinates corresponding to one of the subsets. This is exactly the property that we will consider in the generalisation presented in Section 5.

The multiple-of-8 property

The multiple-of-8 property presented in Lemma 2 is then a direct consequence of the previous theorem. It is derived by combining the theorem with the following proposition, which computes the cardinality of the equivalence classes.

Proposition 1. Let C be an equivalence class with information set K. The cardinality of

C is |C| = 2 |K|-1+d|I|(4-|K|) .
It is always a multiple of 8.

Proof. Since for a given set {p 0 , p 1 } with information set K, we have that ∀k ∈ K, ∀i ∈

I, p 0 i,k = p 1 i,k ,
n = ∆ -1 (D J ) = C∈P 2 (M I +a) ∼ ∆ -1 (D J ) ∩ C .
We know from Theorem 4 that ∆ is constant on the equivalence classes, implying that ∆ -1 (D J ) ∩ C equals either 0 or |C|. In other words, there exists a function δ from

P 2 (M I + a) ∼ into {0, 1} such that ∆ -1 (D J ) ∩ C = δ(C) × |C| . It follows that n = C∈P 2 (M I +a) ∼ δ(C) × |C| ≡ 0 mod 8,
since, by Proposition 1, all equivalence classes have a cardinality divisible by 8.

Influence of the branch number

In [START_REF] Grassi | A new structuraldifferential property of 5-round AES[END_REF], the proof of Lemma 2 uses the fact that the differential branch number of MC denoted by b is maximal and equal to 5 but as we have seen, the branch number b does not have any importance for this lemma. It affects the exact value of n, but not the fact that it is a multiple of 8. Indeed, in the formula given in the proof of Corollary 1 for computing n, we can distinguish between the different equivalent classes according to the size of their information set:

n = C∈P 2 (M I +a) ∼ δ(C) × |C| = 4 h=0 C:|K(C)|=h δ(C) × |C| .
Obviously, for any {p 0 , p 1 } with information set K,

p 0 + p 1 = 3 k=0 i∈I (p 0 i,k + p 1 i,k )MC(e i-k,k ) = k∈K i∈I (p 0 i,k + p 1 i,k )MC(e i-k,k ) ∈ C K , implying,

An alternative proof of Theorem 3

The multiple-of-8 property is a consequence of Theorem 4, which states that the function ∆ is constant over each equivalence class. However, this invariance can be directly used as a distinguishing property. This is actually what is done by the second mixture-differential distinguisher exhibited in [START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF] and detailed in Theorem 3: Theorem 3 is nothing else than the combination of Theorem 4 with the subspace trail given by Lemma 1. Indeed, consider p 0 and p 1 in C i + a with p 0 ≡ (x 0 , x 1 , x 2 , x 3 ) and p 1 ≡ (y 0 , y 1 , y 2 , y 3 ). Then, from Lemma 1, R(p 0 ) and R(p 1 ) belong to the same coset of M i and their decompositions over the basis (MC(e i,0 ), MC(e i-1,1 ), MC(e i-2,2 ), MC(e i-3,3 )) are given by

R(p 0 ) ≡ (S-box(x 0 + a 0,i ), . . . , S-box(x 3 + a 3,i ))
R(p 1 ) ≡ (S-box(y 0 + a 0,i ), . . . , S-box(y 3 + a 3,i )) .

It follows that, if the two pairs {p 0 , p 1 } and {q 0 , q 1 } satisfy one of the relations given in Theorem 3, then {R(p 0 ), R(p 1 )} and {R(q 0 ), R(q 1 )} belong to the same equivalence class for ∼. We then deduce from Theorem 4 that ∆ takes the same value on these two pairs, i.e., R 2 (p 0 ) + R 2 (p 1 ) = R 2 (q 0 ) + R 2 (q 1 ) .

Consequently,

R 2 (p 0 ) + R 2 (p 1 ) ∈ D J if and only if R 2 (q 0 ) + R 2 (q 1 ) ∈ D J .
Moreover, we know from Lemma 1 that

D J R ⇒ C J R ⇒ M J , implying that R 4 (p 0 ) + R 4 (p 1 ) ∈ M J if and only if R 4 (q 0 ) + R 4 (q 1 ) ∈ M J .

Adaptation to a general SPN construction

We provide in this section an extensive version of the equivalence relation and the multipleof-8 property for a more general SPN cipher than the AES. A natural question for this extension is to find out what is the particular property of the subspaces M I for Theorem 4 and Lemma 2 to work and whether these spaces could be replaced by others without altering the result. For a general SPN cipher, we analyse the form of such subspaces with respect to the non-linear layer of the cipher and provide the necessary conditions for their successful combination.

A more general setting for Theorem 4 and Lemma 2

Consider a general SPN cipher working on a state of N words, where the size of each word equals the cipher's S-box size. Suppose that the cipher is iterated for an arbitrary number of rounds and that the round-keys as well as the internal state are represented as word-vectors in K N (and not as matrices). An SPN round R is the composition K • L • S where:

• S is the substitution layer applying an invertible S-box : K → K to each word of the internal state in a certain basis {f 0 , . . . , f N -1 } of K N . It is important to notice that we define S and the basis together.

• L is the linear layer, a bijective F 2 -linear map of K N .

• K is the AddRoundKey operation adding to the internal state a round-key of the same size.

We now want to describe a more general subspace V of K N that could play the role of M I in an adaptation of Theorem 4 to the previously described SPN. As we have already noticed in the previous section, the proof of Theorem 4 relies on the fact that the coordinates of the elements in the input subspace V can be decomposed into several subsets in such a way that each S-box involves only one coordinate subset. This property is captured by the following definition: it requires the existence of a basis of V which can be decomposed over the original basis {f 0 , . . . , f N -1 } by a block-diagonal matrix.

Definition 5. Let V be a subspace of K N . We say that V is compatible with S if there exists a basis of V whose elements written in the basis {f 0 , . . . , f N -1 } form a block-diagonal matrix (with blocks having potentially different dimensions) for a certain order of the elements in both bases. We call such a basis of V a compatibility basis.

Given an arbitrary basis of a subspace V , it is quite easy to check whether V is compatible with S by computing the unique reduced echelon form of the corresponding matrix. If, for a given ordering of the rows, this matrix has a reduced echelon form which is block-diagonal, then V is compatible with S and the reduced echelon form provides a compatibility basis. Otherwise, V is not compatible with S.

We provide now the necessary notation for describing a compatibility basis g of a compatible subspace V . For this notation to be as clear as possible, we first give a representation of g as a collection of column vectors written in the basis f .

                   * • • • * . . . λ 0, ,i . . . 0 0 * • • • * * • • • * 0 . . . λ k, ,i . . . 0 * • • • * * • • • * 0 0 . . . λ h-1, ,i . . . * • • • * 0 0 0                    ← coordinate on f j0+ ← j k + ← j h-1 + ↑ ↑ ↑ g 0,i g k,i g h-1,i i < i 0 i < i k i < i h-1
• We denote by h the number of blocks. The number of basis vectors in the k-th block, 0 ≤ k < h will be denoted by i k . It must obviously hold that

h-1 k=0 i k = dim V .
• The basis of V will be denoted by (g

k,i ) k<h,i<i k ∈ (K N ) dim V , which means that V = vect K (g k,i | k < h, i < i k ).
The index k of a basis element g k,i stands for the block-number, while i represents the position of the vector inside the block k.

• There exist (h + 1) integers j 0 , . . . , j h , with j 0 = 0 and j h ≤ N , such that for all vectors inside the k-th block, all coordinates outside the interval {j k , . . . , j k+1 -1} are zero.

Then, each basis vector g k,i can be written as

g k,i = j k+1 -j k -1 =0 λ k, ,i f j k + (3)
for some λ k, ,i ∈ K with 0 ≤ k < h and 0 ≤ i < i k .

Example 2. For the AES, we have N = 16 and the original basis {f 0 , . . . , f 15 } is formed by the vectors e i,j with i, j ∈ {0, . . . , 3}. Let I = {t 0 , . . . , t r-1 } be a subset of size r of {0, . . . , 3}. Then, M I is compatible with the AES S-box layer. A compatibility basis of M I is given by g k,i = MC(e ti-k,k ) for 0 ≤ k < 4 and 0 ≤ i < r .

Indeed, when the elements of the original basis are ordered by f 4j+i = e i,j for i, j ∈ {0, . . . , 3}, (g 0,0 , . . . , g 3,r-1 ) is obtained by multiplying (f 0 , . . . , f 15 ) by a matrix with h = 4 blocks, all of size r. This comes from the fact that, for j k = 4k, we have

g k,i = MC(e ti-k,k ) = 3 =0 m ,ti-k e ,k = j k+1 -j k -1 =0 m ,ti-k f j k + ,
where m i,j are the coefficients of the 4 × 4 matrix defining MixColumns. This is exactly the block-diagonal form described by (3).

For example, a basis of M 0 can be written as follows, where . corresponds to 0.

                            2 . . . 1 . . . 1 . . . 3 . . . . 1 . . . 1 . . . 3 . . . 2 . . . . 1 . . . 3 . . . 2 . . . 1 . . . . 3 . . . 2 . . . 1 . . . 1                            
From now on, we fix a ∈ K N and a subspace V compatible with S with compatibility basis g. The notion of information set and the equivalence relation between pairs of elements in (a + V ) now involve the decomposition of the elements in V over the basis {g k,i }, where the coordinates corresponding to the same k are gathered together. The next definition adapts the notion of information set to this context. Definition 6 (Information set). Let {p 0 , p 1 } be an unordered pair of elements from V + a, written as

p 0 = h-1 k=0 i k -1 i=0 p 0 i,k g i,k + a and p 1 = h-1 k=0 i k -1 i=0 p 1 i,k g i,k + a
for some (uniquely defined) p 0 i,k , p 1 i,k ∈ K. The information set K of {p 0 , p 1 } is defined as

K = {k ∈ {0, . . . , h -1} | ∃i < i k : p 0 i,k = p 1 i,k } .
Similarly, we define an equivalence relation between pairs of inputs by considering all pairs of elements in (a + V ) obtained by exchanging the sets of coordinates corresponding to the same k.

Definition 7 (Equivalence relation). Let P = {p 0 , p 1 } and Q = {q 0 , q 1 } with p 0 , p 1 , q 0 , q 1 ∈ (V + a). We say that P ∼ Q if:

• {p 0 , p 1 } and {q 0 , q 1 } have the same information set K.

• ∀k ∈ K, ∃b ∈ {0, 1} : ∀i < i k , q 0 i,k = p b i,k and q 1 i,k = p 1-b i,k .
∼ is an equivalence relation on the set of unordered pairs of elements in (V + a).

The next theorem is an adaptation of Theorem 4 to any subspace V compatible with S.

Theorem 5. For any a ∈ K N , the function ∆ operating on unordered pairs of elements in (V + a) and defined by

∆ : {p 0 , p 1 } -→ R(p 0 ) + R(p 1 )
is constant over the equivalence classes for ∼.

Proof. Let P = {p 0 , p 1 }, Q = {q 0 , q 1 } such that P ∼ Q. We have with the previously introduced notation

p 0 = h-1 k=0 i k -1 i=0 p 0 i,k g i,k + a = h-1 k=0 j k+1 -j k -1 =0 i k -1 i=0 p 0 i,k λ i,k, + a j k + f j k +
where the last equality is obtained by replacing each g i,k by its decomposition on the basis (f 0 , . . . , f N -1 ) given by (3). Then

S(p 0 ) = h-1 k=0 j k+1 -j k -1 =0 S-box i k -1 i=0 p 0 i,k λ i,k, + a j k + f j k +
and the difference S(p 0 ) + S(p 1 ) can be written as:

k, S-box i k -1 i=0 p 0 i,k λ i,k, + a j k + + S-box i k -1 i=0 p 1 i,k λ i,k, + a j k + f j k + (4)
It is now clear with Definition 7 and Equation (4) that S(p 0 ) + S(p 1 ) and S(q 0 ) + S(q 1 ) are equal in K N . Therefore,

∆(P ) = R(p 0 ) + R(p 1 ) = K • L • S(p 0 ) + K • L • S(p 1 )
= L(S(p 0 ) + S(p 1 )) by characteristic 2 and linearity of L ; = L(S(q 0 ) + S(q 1 )) by our previous observation ; = ∆(Q).

From this invariance theorem, we can derive, as in the case of the AES distinguisher, some information on the number of unordered pairs {p 0 , p 1 } such that ∆({p 0 , p 1 }) belongs to a given set E. This consequence is similar to the multiple-of-8 property, but the divisibility of this number depends on the structure of the subspace V we consider. This comes from the sizes of the equivalence classes, which are determined in the following proposition.

Proposition 2. Let C be an equivalence class with information set K. The cardinality of

C is |C| = 2 |K|-1+d k ∈K i k .
It is always a multiple of 2 h-1 .

Proof. We have k ∈K (2 d ) i k choices for the shared coordinates in a pair of C. Those coordinates fixed, we have to make for all k ∈ K the choice b = 0 or b = 1, i.e. 2 |K| choices.

Since we are counting unordered pairs, we have 2

|K|-1+d k ∈K i k pairs in C. The exponent |K| -1 + d k ∈K i k is minimal for |K| = h. Indeed,   |K| -1 + d k ∈K i k   = d • dim V -1 - k∈K (d • i k -1)
which obviously decreases as K gets bigger. Hence |C| ≡ 0 mod 2 h-1 .

We then deduce the following generalisation of the multiple-of-8 property.

Corollary 2. Let E be any subset of K N and

n = #{ {p 0 , p 1 } with p 0 , p 1 ∈ (V + a) | R(p 0 ) + R(p 1 ) ∈ E}.
Then n ≡ 0 mod 2 h-1 .

The proof is the same as the proof of Corollary 1, but we detail it for the sake of completeness.

Proof. We denote by P 2 (V + a) the set of all unordered pairs of elements in V + a and by P 2 (M I + a) ∼ the set of all equivalence classes for ∼. Since the equivalence classes form a partition of P 2 (M I + a), we have that

n = ∆ -1 (E) = C∈P 2 (M I +a) ∼ ∆ -1 (E) ∩ C .
We know from Theorem 4 that ∆ is constant on the equivalence classes, implying that

∆ -1 (E) ∩ C = δ(C) × |C| for some function δ from P 2 (M I + a) ∼ into {0, 1}. It follows that n = C∈P 2 (M I +a) ∼ δ(C) × |C| ≡ 0 mod 2 h-1 ,
since, by Proposition 2, all equivalence classes have a cardinality divisible by 2 h-1 .

A new proof of Theorem 2

As a first illustration, we now show that the mixture-differential distinguisher described by Theorem 2 and originally stated in [START_REF] Grassi | Mixture differential cryptanalysis: a new approach to distinguishers and attacks on round-reduced AES[END_REF] can be seen as a direct application of Theorem 5. In this case, the compatible subspace V is vect K (MC(e 0,i ), MC(e 1,i )) = M i ∩ C 0,1 with i ∈ {0, . . . , 3}. The proof of Theorem 2 is then similar to the one presented in Section 4.4, but for a different subspace V . We fix a ∈ M 4 (K) for the rest of this section.

Proof. We define V as the subspace vect K (MC(e 0,i ), MC(e 1,i )). A basis of this subspace is composed of the two column vectors g 0,0 = MC(e 0,i ) and g 1,0 = MC(e 1,i ), whose decomposition over the canonical basis is defined by:

                            2 . 1 . 1 . 3 . . 1 . 1 . 3 . 2 . . . . . . . . . . . . . .
which justifies that V is compatible with the AES substitution layer. Let p 0 , p 1 , q 0 , q 1 ∈ vect K (e 0,i , e 1,i ) + a such that p 0 ≡ (x 0 , x 1 ), p 1 ≡ (y 0 , y 1 ) and q 0 ≡ (x 0 , y 1 ), q 1 ≡ (y 0 , x 1 ) .

Then, R(p 0 ) ≡ (S-box(x 0 +a 0,i ), S-box(x 1 +a 1,i )) and R(p 1 ) ≡ (S-box(y 0 +a 0,i ), S-box(y 1 +a 1,i )) .

This exactly means that {R(p 0 ), R(p 1 )} and {R(q 0 ), R(q 1 )} are equivalent pairs of a coset of V . Theorem 5 then gives that R 2 (p 0 ) + R 2 (p 1 ) = R 2 (q 0 ) + R 2 (q 1 ). Consequently,

R 2 (p 0 ) + R 2 (p 1 ) ∈ D J if and only if R 2 (q 0 ) + R 2 (q 1 ) ∈ D J . Since D J 2R ⇒ M J , this equivalently means that R 4 (p 0 ) + R 4 (p 1 ) ∈ M J if and only if R 4 (q 0 ) + R 4 (q 1 ) ∈ M J .
Obviously the same result also holds for any subspace V formed by the intersection between M i and another subset (unless the matrix defining the corresponding basis has a single block).

Applications

In this section we provide some applications of Theorem 5 to SPN ciphers other than the AES. The goal is to show in practice that the mixture-differential distinguishers and the multiple-of property are not proper to the AES but that they hold for many other SPN constructions. Furthermore, as a result of the adaptation provided in the previous section, the application to other ciphers is almost straightforward. Therefore, we adapt Theorem 1 to the SPN ciphers LED [START_REF] Guo | The LED block cipher[END_REF], Midori [BBI + 15], Klein [START_REF] Gong | KLEIN: A new family of lightweight block ciphers[END_REF] and Skinny [BJK + 16] and discuss why the result does not adapt to Crypton [START_REF] Hoon | A revised version of Crypton -Crypton v1.0[END_REF] or Prince [BCG + 12].

In practice, for finding a multiple-of property for some cipher, two conditions must be met. The first one is the existence of a subspace V compatible with the substitution layer for the multiple-of property through one round to hold. This condition is described through Corollary 2, which can be seen as a general replacement of Lemma 2. The second condition is the existence of exact subspace trails covering some rounds before and after the central round on which the multiple-of property is found. For this, we need a result equivalent to Lemma 1 for each analysed cipher. For searching for subspace trails that can replace the subspaces D I , C I and M I in Lemma 1, we use Algorithm 1 proposed by Leander, Tezcan and Wiemer in [START_REF] Leander | Searching for subspace trails and truncated differentials[END_REF].

Midori, Klein, Skinny and Crypton all follow our general SPN description and operate on a state composed of N = 16 words, where the size of each word (i.e. the dimension of the S-box alphabet) is d ∈ {4, 8}. LED and a part of Prince also follow this general description with d = 4. We study for each cipher the exact subspace trails given by Algorithm 1 in [START_REF] Leander | Searching for subspace trails and truncated differentials[END_REF] and the compatibility of the last subspace of the trail with the substitution layer in the sense of Definition 5.

The cases of AES, LED, Midori, Klein and Skinny

AES. First of all, a natural idea would be to use our adaptation with a longer subspace trail than the ones of Lemma 1 for the AES. However, as shown in [START_REF] Leander | Searching for subspace trails and truncated differentials[END_REF]

    0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0    
to the columns of the state. The branch number of the above matrix is 4. SC is the so-called ShuffleCell permutation that permutes the words of the state as follows:

(s 0 , s 1 , . . . , s 15 ) ← (s 0 , s 7 , s 14 , s 9 , s 5 , s 2 , s 11 , s 12 , s 15 , s 8 , s 1 , s 6 , s 10 , s 13 , s 4 , s 3 ) where the words are numbered column-wise. As in Section 2, we define the following subspaces depending on the linear components of Midori's round-function. If i ∈ {0, 1, 2, 3},

C i = vect K (e 0,i , e 1,i , e 2,i , e 3,i ), D Mi i = SC -1 (C i ), M Mi i = L Mi (C i ).
Applying Algorithm 1 in [START_REF] Leander | Searching for subspace trails and truncated differentials[END_REF] gives that the longest exact subspace trails are the two-round trails of the form

D Mi I R Mi ⇒ C I R Mi ⇒ M Mi I .
Besides, M Mi I has a basis whose matrix is block-diagonal with four blocks. For example, M Mi 0 has a basis whose representation as a collection of column vectors is given in Figure 1a. We then have by Corollary 2 that for all I, J ⊆ {0, 1, 2, 3}, a ∈ (F by (f i ) 0≤i<16 and then define the following subspaces for i ∈ {0, 1} :

2 d ) 16 , #{{p 0 , p 1 } with p 0 , p 1 ∈ D Mi I + a | R 5 Mi (p 0 ) + R 5 Mi (p 1 ) ∈ M Mi J } ≡ 0 mod
C i = vect K (f k | 0 ≤ k < 8), D Kl i = RN -1 (C i ), M Kl i = L Kl (C i ).
Algorithm 1 in [START_REF] Leander | Searching for subspace trails and truncated differentials[END_REF] gives that the longest exact subspace trails are two-round trails of the form Besides, M Kl i has a basis whose matrix is block-diagonal with two blocks. For example, M Kl 0 has a basis whose representation as a collection of column vectors is given in Figure 1b. We then have by Corollary 2 for all i, j ∈ {0, 1}, a ∈ (F 2 4 ) 16 , #{{p 0 , p1 } with p 0 , p 1 ∈ D

D Kl i R Kl ⇒ C i R Kl ⇒ M Kl i .                             0
    1 0 1 1 1 0 0 0 0 1 1 0 1 0 1 0    
Applying here again Algorithm 1 of [START_REF] Leander | Searching for subspace trails and truncated differentials[END_REF] gives as result that the longest exact subspace trails are two-round-long. There are 1294 such trails. A Gauss elimination on the last subspace of each trail gives that among these trails, 1282 end with a subspace compatible with the substitution layer. This allows to adapt Theorem 1 for 5-round Skinny, concluding that 5-round Skinny always has the multiple-of-2 h-1 property. However, it is interesting to note here that depending on the trail, the value of h varies. More precisely, h ∈ {2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14}. We now give an example of such a distinguisher for h = 3. For Skinny, the internal state is classically represented as an element of M 4 (K). We first exhibit two 2-round subspace trails, U i

R Sk ⇒ V i R Sk
⇒ W i for i ∈ {0, 1} where U 0 = vect K (e 1,1 , e 1,2 , e 1,3 , e 3,1 , e 3,3 ), V 0 = L Sk (U 0 ), W 0 = L Sk (V 0 ) and U 1 = vect K (e 0,3 , e 1,0 , e 1,2 , e 1,3 , e 2,1 , e 2,3 , e 3,0 , e 3,1 , e 3,2 , e 3,3 ), V 1 = L Sk (U 1 ), W 1 = L Sk (V 1 ) .

Moreover, W 0 is compatible with h = 3. One of its compatibility basis is represented in Figure 1c. We then have that #{{p 0 , p 1 } with p 0 , p 1 ∈ U 0 + a | R 5

Sk (p 0 ) + R 5 Sk (p 1 ) ∈ W 1 } ≡ 0 mod 4 .

Finally, we recall that this section only aims at giving examples of how Theorem 1 can be adapted to other SPN ciphers and does not claim new cryptanalytic results. In particular, this property on 5 rounds does not threaten the overall security of the cipher that is composed of 32 rounds.

The cases of Crypton and Prince

Crypton The block cipher Crypton [START_REF] Hoon | CRYPTON: a New 128-bit Block Cipher[END_REF], designed by Lim in 1998 was among the candidates to the NIST AES competition. It has a structure very similar to the one of AES and this is why we considered it as a natural candidate for the multiple-of property. Indeed, the round function of Crypton is naturally decomposed as R Cr = L Cr • S Cr , where L Cr is the composition of a byte transposition of columns into rows with respect to the anti-diagonal of the internal state and a permutation at the bit level applied column-wise. Algorithm 1 in [START_REF] Leander | Searching for subspace trails and truncated differentials[END_REF] gives two-round exact subspace trails. However, the problem here is that L Cr is only F 2 -linear and not F 2 8 -linear and this implies that the last subspaces in the subspace trails are not F 2 8 -vector subspaces and cannot verify the compatibility hypothesis of Corollary 2. As a consequence, Theorem 1 cannot be adapted to 5-round Crypton. However, Crypton has 1-round exact subspace trails that end with subspaces of the form vect K (e i,j |i ∈ I, j ∈ {0, 1, 2, 3}). Those subspaces are obviously compatible with the substitution layer, and we have a 4-round version of Theorem 1 for Crypton: the first round is a 1-round exact subspace trail ending with a compatible subspace, the second round exploits Corollary 2 and the last two rounds are a 2-round exact subspace trail. Indeed, the trail used after Corollary 2 does not need to end with a compatible subspace.

Prince Prince [BCG + 12] is a block cipher proposed by Borghoff et al. in 2012 with a specific structure named α-reflection. The main parts of this structure follow the SPN construction. As for Crypton, the round function is defined as R Pr = L Pr • S Pr with L Pr a F 2 -linear map which is not F 2 4 -linear. Again, as for Crypton, Prince exhibits two-round exact subspace trails but the last subspaces of those trails are not F 2 4 -linear subspaces. We then cannot mount a 5-round distinguisher by adapting Theorem 1. However, one-round exact subspace trails ending with compatible diagonal F 2 4 -linear subspaces allow to have a 4-round version of Theorem 1 for Prince.

Conclusion

We have presented a general result which allows cryptanalysts to search for mixturedifferential distinguishers, or multiple-of properties, in a systematic way, for any SPN. This result then avoids the redundant proofs which were previously needed for each new occurrence of these distinguishing properties. Also, it highlights the properties of the ciphers which have to be taken into account for establishing the existence of such distinguishers and it shows that mixture-differential distinguishers directly apply to a more general class of SPN than what was previously believed. As shown in the previous examples, all these distinguishing properties can be exhibited by combining our framework with the search for subspace trails, which is investigated in [START_REF] Leander | Searching for subspace trails and truncated differentials[END_REF]. Since our result, exploiting an appropriate equivalence relation, applies in many situations, it appears that the main limitation for finding efficient distinguishers is the existence of long subspace trails for the ciphers.

  we have (2 d ) |I|×(4-|K|) choices for the shared coordinates in a pair of C. Those coordinates fixed, we have to make for all k ∈ K the choice b = 0 or b = 1, i.e. 2 |K| choices. Since we are counting unordered pairs, we have 2 |K|-1+d|I|(4-|K|) elements in C. Obviously, the exponent |K| -1 + d|I|(4 -|K|) is minimal for |K| = 4. We deduce that |K| -1 + d|I|(4 -|K|) ≥ 3 , leading to |C| ≡ 0 mod 8. By combining Proposition 1 and Theorem 4, we deduce the following corollary which generalises Lemma 2 in the sense that we are not restricted to the case |I| = 1 as in Lemma 2. Corollary 1. Let n = #{ {p 0 , p 1 } with p 0 , p 1 ∈ M I + a | R(p 0 ) + R(p 1 ) ∈ D J }. Then, n ≡ 0 mod 8. Proof. Let P 2 (M I + a) denote the set of all unordered pairs of elements in M I + a, and P 2 (M I + a) ∼ denote the set of all equivalence classes for ∼. Since the equivalence classes form a partition of P 2 (M I + a), we have that

  by Lemma 1, that (R(p 0 ) + R(p 1 )) belongs to M K . However, it has been proved in [GRR16, Lemma 5] that, for any two sets I, J ⊆ {0, . . . , 3} such that |I| + |J| < b where b is the branch number of MC, we have D I ∩ M J = {0}. It then follows that, if K has size h < b -|J|, then R(p 0 ) + R(p 1 ) ∈ D J unless p 0 = p 1 . We can then express the influence of the branch number on n with the formula n = 4 h=b-|J| C:|K(C)|=h δ(C) × |C| , which does not affect the multiple-of-8 property.

  , the exact trails D I R ⇒ C I R ⇒ M I are the longest possible trails for this cipher, which means that Corollary 2 cannot give any improvement for the AES. LED. LED is a lightweight block cipher proposed by Guo et al. [GPPR11]. Its round function has the same structure as the AES and exhibits the same two-round subspace trails D I R ⇒ C I R ⇒ M I , leading to a 5-round distinguisher. Midori. Midori is a lightweight block cipher designed by Banik et al. [BBI + 15], optimized with respect to the energy consumption. The round function R Mi is the composition L Mi • S Mi , where L Mi = MC Mi • SC. MC Mi is the MixColumns operation, applying the binary involutive matrix

Figure 2 :

 2 Figure 2: Compatibility basis representations for Midori (left) and Klein (right)

  Klein is a lightweight block cipher proposed in 2011 by Gong et al.[START_REF] Gong | KLEIN: A new family of lightweight block ciphers[END_REF]. The round function of Klein R Kl = L Kl • S Kl can be seen as the application of a non-linear layer S Kl followed by a linear layer L Kl . The linear layer L Kl = MN • RN is the composition of RN, standing for the RotateNibbles permutation rotating the state two bytes to the left and MN, standing for the MixNibbles permutation. This last operation applies the AES MixColumns transformation to each half of the state. We denote the canonical basis of F 16 2 4

8, which gives a 5-round distinguisher similar to the AES-one. It is worth noticing that the property holds even if the branch number of the MixColumns operation in Midori is only 4.

Klein.

  Skinny is a family of tweakable lightweight block ciphers, designed in 2016 by Beierle et al. [BJK + 16]. The round function follows a classical SPN construction R Sk = L Sk • S Sk where L Sk = MC Sk • SR Sk . The operation SR Sk is similar to the AES ShiftRows operation, with the only difference that the shift is performed to the right. MC Sk is the MixColumns operation, where each column of the state is multiplied by the following binary matrix of branch number 2:

Kl i + a | R 5 Kl (p 0 ) + R 5 Kl (p 1 ) ∈ M Kl j } ≡ 0 mod 2 ,

meaning that Klein has a multiple-of-2 property for 5 rounds. Note that even if we get only a multiple of 2 in the case of Klein, as the pairs are not ordered, this can still be considered as a distinguishing property.

Skinny.

In

2011, biclique attacks were applied against all full-round versions of the AES[START_REF] Bogdanov | Biclique cryptanalysis of the full AES[END_REF]. This technique permitted to reduce the exhaustive key-search by a few bits. However, we consider here that this kind of attack is a form of accelerated key-search.

. . . . . 1 . . . 1 1 . . . 1 . . . . . . 1 . . . . 1 . . . . . . . . . 1 . . . . . 1 . . . . . 1 . . . 1 1 . . . 1 . . . . . . . . . . .

Acknowledgments

We would like to thank Hadi Soleimany and the anonymous reviewers for their helpful comments and suggestions.


V Shutyaev
email: victor.shutyaev@mail.ru

F.-X Le Dimet

E Parmuzin

Sensitivity of response functions in variational data assimilation for joint parameter and initial state estimation

Keywords: Data assimilation, optimal control, parameter estimation, sensitivity, sea thermodynamics model

The problem of variational data assimilation for a nonlinear evolution model is formulated as an optimal control problem to find simultaneously unknown parameters and initial state of the model. A response function is considered as a functional of the optimal solution after assimilation. The sensitivity of the response function to the observation data is studied. The gradient of the response function with respect to observations is related to the solution of a non-standard problem involving the coupled system of direct and adjoint equations. Based on the Hessian of the original cost function, the solvability of the non-standard problem is studied. An algorithm to compute the gradient of the response function with respect to observation data is formulated and justified. A numerical example is presented for variational data assimilation problem for the Baltic Sea thermodynamics model.

Introduction

The methods of variational data assimilation have become a very important tool for state observation and parameter estimation for geophysical models. The problems of variational data assimilation can be formulated as optimal control problems (e.g. [START_REF] Lions | Contrôle Optimal des Systèmes Gouvernés par des Équations aux Dérivées Partielles[END_REF]- [START_REF] Carrassi | Data assimilation in the geosciences: An overview of methods, issues, and perspectives[END_REF]) to find unknown model parameters such as initial and boundary conditions, right-hand sides in the model equations, distributed coefficients, based on minimization of the cost function related to observations. A necessary optimality condition reduces an optimal control problem to an optimality system which involves the model equations, the adjoint problem, and input data functions. The optimal solution depends on the observation data, which may contain uncertainties, and for the forecasts it is very important to study the sensitivity of the optimal solution and its functionals with respect to observation errors [START_REF] Baker | Observation and background adjoint sensitivity in the adaptive observation-targeting problem[END_REF].

The necessary optimality condition is related to the gradient of the original cost function, thus to study the sensitivity of the optimal solution, one should differentiate the optimality system with respect to observations. In this case, we come to the so-called second-order adjoint problem [START_REF] Dimet | Second-order information in data assimilation[END_REF]. The first studies of sensitivity of the response functions after assimilation with the use of second-order adjoint were done by [START_REF] Dimet | Sensitivity analysis in variational data assimilation[END_REF] for variational data assimilation problem aimed at restoration of initial condition, where sensitivity with respect to model parameters was considered. The equations of the forecast sensitivity to observations in a four-dimensional (4D-Var) data assimilation were derived by [START_REF] Daescu | On the sensitivity equations of four-dimensional variational (4D-Var) data assimilation[END_REF]. Based on these results, a practical computational approach was given by [START_REF] Cioaca | Efficient methods for computing observation impact in 4D-Var data assimilation[END_REF] to quantify the effect of observations in 4D-Var data assimilation.

Sensitivity of the optimal solution is related to its statistical properties (see [START_REF] Gejadze | On analysis error covariances in variational data assimilation[END_REF]- [START_REF] Shutyaev | Optimal solution error covariance in highly nonlinear problems of variational data assimilation[END_REF]). General sensitivity analysis in variational data assimilation with respect to observations for a nonlinear dynamic model was given by [START_REF] Shutyaev | Sensitivity with respect to observations in variational data assimilation[END_REF] to control the initial-value function. The paper [START_REF] Shutyaev | Sensitivity analysis with respect to observations in variational data assimilation for parameter estimation[END_REF] presented the sensitivity analysis with respect to observations in variational data assimilation aimed at restoration of unknown parameters of a dynamic model.

This paper generalizes the results of [START_REF] Shutyaev | Sensitivity with respect to observations in variational data assimilation[END_REF] and [START_REF] Shutyaev | Sensitivity analysis with respect to observations in variational data assimilation for parameter estimation[END_REF] and presents the sensitivity analysis with respect to observations in variational data assimilation aimed at simultaneous restoration of unknown parameters and initial state of a dynamic model. The problems of parameter estimation are common inverse problems considered in geophysics and in engineering applications (see [START_REF] Alifanov | Extreme Methods for Solving Ill-posed Problems with Applications to Inverse Heat Transfer Problems[END_REF]- [START_REF] Agoshkov | The study and numerical solution of some inverse problems in simulation of hydrophysical fields in water areas with liquid boundaries[END_REF]). Last years an interest is rising to the joint initial state and parameter estimation using 4D-Var ( [START_REF] Dee | Bias and data assimilation[END_REF] - [START_REF] Herzog | Optimal sensor placement for joint parameter and state estimation problems in large-scale dynamical systems with applications to thermo-mechanics[END_REF]).

We consider a dynamic formulation of variational data assimilation problem for joint parameter and initial state estimation in a continuous form, but the presented sensitivity analysis formulas with respect to observations do not follow from our previous results [START_REF] Shutyaev | Sensitivity with respect to observations in variational data assimilation[END_REF] and [START_REF] Shutyaev | Sensitivity analysis with respect to observations in variational data assimilation for parameter estimation[END_REF] and constitute a novelty of this paper. The main contribution of the paper, as compared to [START_REF] Shutyaev | Sensitivity with respect to observations in variational data assimilation[END_REF] and [START_REF] Shutyaev | Sensitivity analysis with respect to observations in variational data assimilation for parameter estimation[END_REF], is a derivation of new formulas for the gradient of a response function with respect to observations in variational data assimilation problem aimed at joint parameter and initial state estimation for a general nonlinear dynamic model. This paper is organized as follows. In section 2, we give the statement of the variational data assimilation problem for a nonlinear evolution model to estimate simultaneously the model parameters and the initial state. In Section 3, sensitivity of the response function after assimilation with respect to observations is studied, and the theorem is proved to relate the gradient to the solution of a non-standard problem. An algorithm to compute the gradient of the response function is formulated, based on an operator equation involving the Hessian of the original cost function, and the the solvability of the non-standard problem is studied. In Section 4, we consider a simple example with known manufactured exact solution and present some numerical results.

Section 5 presents an application of the theory to the data assimilation problem for the Baltic Sea thermodynamics model with a numerical example. The main results are discussed in the Conclusions.

Statement of the problem

Consider the mathematical model of a physical process that is described by the nonlinear evolution problem

∂ϕ ∂t = F (ϕ, λ) + f, t ∈ (0, T ) ϕ t=0 = u, (1)
where the initial state u is supposed to be from a Hilbert space X, the unknown function ϕ = ϕ(t) belongs to Y = L 2 (0, T ; X) with the norm ϕ Y = (ϕ, ϕ)

1/2 Y = ( T 0 ϕ(t) 2 X dt) 1/2 , F is a nonlinear operator mapping Y × Y p into Y , Y p is a Hilbert space (space of model parameters),
f ∈ Y . We suppose that for given u ∈ X, f ∈ Y and λ ∈ Y p there exists a unique solution ϕ ∈ Y to [START_REF] Lions | Contrôle Optimal des Systèmes Gouvernés par des Équations aux Dérivées Partielles[END_REF] with ∂ϕ ∂t ∈ Y . The function λ is an unknown model parameter, and we suppose that the initial state u is also unknown.

We introduce the cost function as a functional on X × Y p in the form

J(u, λ) = 1 2 (V 1 (u-u b ), u-u b ) X + 1 2 (V 2 (λ-λ b ), λ-λ b ) Yp + 1 2 (V 3 (Cϕ-ϕ obs ), Cϕ-ϕ obs ) Y obs , (2) where u b ∈ X, λ b ∈ Y p are prior (background) functions, ϕ obs ∈ Y obs is a prescribed function (observational data), Y obs is a Hilbert space (observation space), C : Y → Y obs is a linear bounded operator (observation operator), V 1 : X → X, V 2 : Y p → Y p and V 3 : Y obs → Y obs are symmetric positive definite bounded operators. Usually, V 1 , V 2 , V 3 are
chosen as inverse covariance operators of observation and background errors [START_REF] Baker | Observation and background adjoint sensitivity in the adaptive observation-targeting problem[END_REF][START_REF] Smith | Data assimilation for state and parameter estimation: application to morphodynamic modelling[END_REF].

Let us consider the following data assimilation problem with the aim to find the initial value u and the parameter λ:

for given f ∈ Y, ϕ obs ∈ Y obs , u b ∈ X, λ b ∈ Y p , find u ∈ X, λ ∈ Y p and
ϕ ∈ Y such that they satisfy (1), and on the set of solutions to (1), the functional J(u, λ) takes the minimum value, i.e.

       ∂ϕ ∂t = F (ϕ, λ) + f, t ∈ (0, T ) ϕ t=0 = u, J(u, λ) = inf w∈X,v∈Yp J(w, v). (3)
We suppose that the solution of (3) exists. Let us note that the solvability of the parameter estimation problems (or identifiability) has been addressed, e.g., in [START_REF] Chavent | Local stability of the output least square parameter estimation technique[END_REF], [START_REF] Navon | Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography[END_REF]. To derive the optimality system, we assume the solution ϕ and the operator F (ϕ, λ) in (1)-(2) are regular enough, and for w ∈ X, v ∈ Y p find the gradient of the functional J with respect to u and λ:

J ′ u (u, λ)w = (V 1 (u -u b ), w) X + (C * V 3 (Cϕ -ϕ obs ), φ) Y , (4)
J ′ λ (u, λ)v = (V 2 (λ -λ b ), v) Yp + (C * V 3 (Cϕ -ϕ obs ), φ) Y , (5)
where φ is the solution to the problem:

∂φ ∂t = F ′ ϕ (ϕ, λ)φ + F ′ λ (ϕ, λ)v, t ∈ (0, T ), φ t=0 = 0, (6)
and φ is the solution to the problem:

∂ φ ∂t = F ′ ϕ (ϕ, λ) φ, t ∈ (0, T ), φ t=0 = w. (7)
Here [START_REF] Marchuk | Adjoint Equations and Perturbation Algorithms in Nonlinear Problems[END_REF] with respect to ϕ and λ, correspondingly, and C * is the adjoint operator to C defined by (Cϕ, ψ)

F ′ ϕ (ϕ, λ) : Y → Y, F ′ λ (ϕ, λ) : Y p → Y are the Fréchet derivatives of F [
Y obs = (ϕ, C * ψ) Y , ϕ ∈ Y, ψ ∈ Y obs .
Let us introduce the adjoint operator (F ′ ϕ (ϕ, λ)) * : Y → Y and consider the adjoint problem:

∂ϕ * ∂t + (F ′ ϕ (ϕ, λ)) * ϕ * = C * V 3 (Cϕ -ϕ obs ), t ∈ (0, T ) ϕ * t=T = 0. (8)
The problem ( 8) is adjoint with respect to the linearized (tangent linear) problems ( 6), [START_REF] Baker | Observation and background adjoint sensitivity in the adaptive observation-targeting problem[END_REF], therefore, it is also linear in ϕ * , however, it is still nonlinear in ϕ.

In the below consideration, we assume that the direct and adjoint linear problems of the form

∂φ ∂t -F ′ ϕ (ϕ, λ)φ = p, t ∈ (0, T ) φ t=0 = q, - ∂φ * ∂t -(F ′ ϕ (ϕ, λ)) * φ * = g, t ∈ (0, T ) φ * t=T = 0
with p, g ∈ Y, q ∈ X have the unique solutions φ, φ * ∈ Y and ∂φ ∂t , ∂φ * ∂t ∈ Y . From ( 4)-( 8) we get

J ′ u (u, λ)w = (V 1 (u -u b ), w) X -(ϕ * t=0 , w) X , (9)
J ′ λ (u, λ)v = (V 2 (λ-λ b ), v) Yp -(ϕ * , F ′ λ (ϕ, λ)v) Y = (V 2 (λ-λ b ), v) Yp -((F ′ λ (ϕ, λ)) * ϕ * , v) Yp , (10) where
(F ′ λ (ϕ, λ)) * : Y → Y p is the operator adjoint to F ′ λ (ϕ, λ).
Thus, the gradient of J is defined by

J ′ u (u, λ) = V 1 (u -u b ) -ϕ * t=0 , J ′ λ (u, λ) = V 2 (λ -λ b ) -(F ′ λ (ϕ, λ)) * ϕ * .
The necessary optimality condition [START_REF] Lions | Contrôle Optimal des Systèmes Gouvernés par des Équations aux Dérivées Partielles[END_REF] is gradJ = 0, therefore, J ′ u (u, λ) = 0, J ′ u (u, λ) = 0. From (3)- [START_REF] Daescu | On the sensitivity equations of four-dimensional variational (4D-Var) data assimilation[END_REF] we obtain the optimality system:

∂ϕ ∂t = F (ϕ, λ) + f, t ∈ (0, T ), ϕ t=0 = u, (11)
∂ϕ * ∂t + (F ′ ϕ (ϕ, λ)) * ϕ * = C * V 3 (Cϕ -ϕ obs ), t ∈ (0, T ) ϕ * t=T = 0, (12)
V 1 (u -u b ) -ϕ * t=0 = 0, (13)
V 2 (λ -λ b ) -(F ′ λ (ϕ, λ)) * ϕ * = 0. (14)
We suppose that the system (11)-( 14) has a unique solution ϕ, ϕ * ∈ Y, u ∈ X, λ ∈ Y p . The system (11)-( 14) may be considered as a generalized model of the form A(U ) = 0 with the state variable U = (ϕ, ϕ * , u, λ), and it contains the information on the observation data ϕ obs ∈ Y obs . Below we study the sensitivity of functionals of the optimal solution with respect to the observation data.

Sensitivity of response functions with respect to observations

In many applications the observation data cannot be measured precisely, and therefore, it is important to be able to estimate the impact of uncertainties in observations on the outputs of the model after assimilation. Such outputs may be response functions considered as functionals of the optimal solution.

We introduce a response function G(ϕ, u, λ), which is supposed to be a real-valued function and can be considered as a functional on Z = Y × X × Y p . We are interested in the sensitivity of G with respect to ϕ obs , with ϕ, u and λ obtained from the optimality system ( 11)- [START_REF] Gejadze | On gauss-verifiability of optimal solutions in variational data assimilation problems with nonlinear dynamics[END_REF]. By definition, the sensitivity is defined by the gradient of G with respect to ϕ obs :

dG dϕ obs = ∂G ∂ϕ ∂ϕ ∂ϕ obs + ∂G ∂λ ∂λ ∂ϕ obs + ∂G ∂u ∂u ∂ϕ obs , (15)
where Let δϕ obs be a perturbation on ϕ obs , then we obtain from the optimality system ( 11)-( 14):

∂G ∂ϕ : Z → Y, ∂G ∂λ : Z → Y p , ∂G ∂u : Z → X,
∂δϕ ∂t = F ′ ϕ (ϕ, λ)δϕ + F ′ λ (ϕ, λ)δλ, t ∈ (0, T ) δϕ t=0 = δu, (16)
     - ∂δϕ * ∂t -(F ′ ϕ (ϕ, λ)) * δϕ * -(F ′′ ϕϕ (ϕ, λ)δϕ) * ϕ * = (F ′′ ϕλ (ϕ, λ)δλ) * ϕ * -C * V 3 (Cδϕ -δϕ obs ), δϕ * t=T = 0, (17)
V 1 δu -δϕ * t=0 = 0, (18)
V 2 δλ -(F ′′ λϕ (ϕ, λ)δϕ) * ϕ * -(F ′′ λλ (ϕ, λ)δλ) * ϕ * -(F ′ λ (ϕ, λ)) * δϕ * = 0, (19)
and dG dϕ obs , δϕ obs

Y obs = ∂G ∂ϕ , δϕ Y + ∂G ∂λ , δλ Yp + ∂G ∂u , δu X , (20)
where δϕ, δϕ * , δλ, δu are the solutions of ( 16)- [START_REF] Alifanov | Extreme Methods for Solving Ill-posed Problems with Applications to Inverse Heat Transfer Problems[END_REF].

The following statement is valid. Theorem 1. Let P 1 , P 2 ∈ Y, P 3 ∈ Y p , P 4 ∈ X be the solutions of the following system of equations

-∂P 1 ∂t -(F ′ ϕ (ϕ, λ)) * P 1 -(F ′′ ϕϕ (ϕ, λ)P 2 ) * ϕ * = (F ′′ λϕ (ϕ, λ)P 3 ) * ϕ * -C * V 3 CP 2 + ∂G ∂ϕ , P 1 t=T = 0, (21)
∂P 2 ∂t -F ′ ϕ (ϕ, λ)P 2 -F ′ λ (ϕ, λ)P 3 = 0, t ∈ (0, T ) P 2 t=0 -P 4 = 0, (22)
V 1 P 4 -P 1 t=0 = ∂G ∂u , (23)
V 2 P 3 -(F ′′ ϕλ (ϕ, λ)P 2 ) * ϕ * -(F ′′ λλ (ϕ, λ)P 3 ) * ϕ * -(F ′ λ (ϕ, λ)) * P 1 = ∂G ∂λ , (24)
where ϕ, ϕ * ∈ Y, u ∈ X, λ ∈ Y p are the solution of the optimality system ( 11)- [START_REF] Gejadze | On gauss-verifiability of optimal solutions in variational data assimilation problems with nonlinear dynamics[END_REF]. Then the gradient of G with respect to ϕ obs is given by

dG dϕ obs = V 3 CP 2 . ( 25
)
Proof of this Theorem is presented in the Appendix. We obtain a coupled system of two differential equations ( 21) and ( 22) of the first order with respect to time, with additional conditions ( 23)- [START_REF] Bocquet | Parameter-field estimation for atmospheric dispersion: application to the Chernobyl accident using 4D-Var, Q[END_REF]. To study this non-standard problem ( 21)-( 24) with mutually dependent initial conditions for P 1 , P 2 , we reduce it to a single operator equation involving the Hessian of the original cost function.

Let us introduce the auxiliary variables v = P 3 ∈ Y p , w = P 4 ∈ X and rewrite the nonstandard problem ( 21)-( 24) in an equivalent form:

∂P 2 ∂t -F ′ ϕ (ϕ, λ)P 2 = F ′ λ (ϕ, λ)v, t ∈ (0, T ) P 2 t=0 = w, (26)
-∂P 1 ∂t -(F ′ ϕ (ϕ, λ)) * P 1 -(F ′′ ϕϕ (ϕ, λ)P 2 ) * ϕ * = (F ′′ λϕ (ϕ, λ)v) * ϕ * -C * V 3 CP 2 + ∂G ∂ϕ , P 1 t=T = 0, (27)
V 1 w -P 1 t=0 = ∂G ∂u , (28)
V 2 v -(F ′′ ϕλ (ϕ, λ)P 2 ) * ϕ * -(F ′′ λλ (ϕ, λ)v) * ϕ * -(F ′ λ (ϕ, λ)) * P 1 = ∂G ∂λ (29)
with the four unknowns: w ∈ X, v ∈ Y p , P 1 , P 2 ∈ Y . Let us write ( 26)-( 29) in the form of an operator equation for U = (w, v) T . We define the operator H : X × Y p → X × Y p , which acts on U belonging to X × Y p , by the successive solution of the following problems:

∂φ ∂t -F ′ ϕ (ϕ, λ)φ = F ′ λ (ϕ, λ)v, t ∈ (0, T ) φ t=0 = w, (30)
- ∂φ * ∂t -(F ′ ϕ (ϕ, λ)) * φ * -(F ′′ ϕϕ (ϕ, λ)φ) * ϕ * = (F ′′ λϕ (ϕ, λ)w) * ϕ * -C * V 3 Cφ, φ * t=T = 0, (31)
HU = V 1 w-φ * t=0 , V 2 v-(F ′′ ϕλ (ϕ, λ)φ) * ϕ * -(F ′′ λλ (ϕ, λ)w) * ϕ * -(F ′ λ (ϕ, λ)) * φ * T , (32)
where λ, u, ϕ and ϕ * are the solutions of the optimality system ( 11)- [START_REF] Gejadze | On gauss-verifiability of optimal solutions in variational data assimilation problems with nonlinear dynamics[END_REF]. It is easily seen that ( 26)-( 29) is equivalent to the following equation in X × Y p :

HU = F (33)
with F ∈ X × Y p defined by

F = ∂G ∂u + φ * t=0 , ∂G ∂λ + (F ′ λ (ϕ, λ)) * φ * T , (34)
where φ * ∈ Y is the solution to the adjoint problem:

   - ∂ φ * ∂t -(F ′ ϕ (ϕ, λ)) * φ * = ∂G ∂ϕ , t ∈ (0, T ) φ * t=T = 0. ( 35
)
It is easy to make sure that the operator H defined by ( 30)-( 32) is the Hessian of the original functional J considered on the optimal solution u, λ of the problem ( 11)-( 14): J ′′ (u, λ) = H.

Lemma 1. Under the assumption that H is positive definite, the operator equation ( 33) is well posed: for every F ∈ X × Y p there exists a unique solution U ∈ X × Y p and the estimate is valid:

U X×Yp ≤ c F X×Yp , c = const > 0. ( 36
)
Proof. If the operator H is positive definite, then for any

U ∈ X × Y p (HU, U ) X×Yp ≥ γ(U, U ) X×Yp , γ = const > 0.
Hence,

HU X×Yp ≥ γ U X×Yp , (37)
and it means that equation ( 33) is uniquely and correctly solvable in X × Y p [START_REF] Vainberg | Variational Methods for the Study of Nonlinear Operators[END_REF].

By definition, H is self-adjoint, i.e. H * = H. Then, the adjoint equation is also correctly solvable, which implies that equation ( 33) is everywhere solvable [START_REF] Vainberg | Variational Methods for the Study of Nonlinear Operators[END_REF], i.e. for every F ∈ X ×Y p there exists a unique solution U ∈ X × Y p .

Let U ∈ X × Y p be the solution of [START_REF] Navon | Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography[END_REF] with the right-hand side F , then [START_REF] Marchuk | Adjoint Equations and Analysis of Complex Systems[END_REF] gives ( 36) with c = 1/γ. The lemma is proved.

Therefore, under the assumption that J ′′ (u, λ) is positive definite on the optimal solution, the non-standard problem ( 21)-( 24) has a unique solution

P 1 , P 2 ∈ Y, P 3 ∈ Y p , P 4 ∈ X.
From the above consideration, we come to the following algorithm to compute the gradient of the response function G:

1) For ∂G ∂λ ∈ Y p , ∂G ∂ϕ ∈ Y, ∂G ∂u ∈ X solve the adjoint problem    - ∂ φ * ∂t -(F ′ ϕ (ϕ, λ)) * φ * = ∂G ∂ϕ , t ∈ (0, T ) φ * t=T = 0 (38)
and put

F = ∂G ∂u + φ * t=0 , ∂G ∂λ + (F ′ λ (ϕ, λ)) * φ * T .
2) Find U = (w, v) T by solving HU = F with the Hessian of the original functional J defined by ( 30)- [START_REF] Chavent | Local stability of the output least square parameter estimation technique[END_REF].

3) Solve the direct problem

∂P 2 ∂t -F ′ ϕ (ϕ, λ)P 2 = F ′ λ (ϕ, λ)v, t ∈ (0, T ) P 2 t=0 = w. ( 39
)
4) Compute the gradient of the response function as

dG dϕ obs = V 3 CP 2 . ( 40
)
The last formula allows us to estimate the sensitivity of the response functions related to the optimal solution after assimilation, with respect to observation data.

Simple example

Let us consider a simple evolution problem for the ordinary differential equation

   dϕ dt + aϕ = λg, t ∈ (0, T ) ϕ t=0 = u, (41)
where u ∈ R; a, λ ∈ R, g = g(t) ≥ 0. Here, in the notations of section 2, we have

X = R, Y = L 2 (0, T ), Y p = R, F (ϕ, λ) = -aϕ + λg, f = 0.
Let us formulate the data assimilation problem to find the initial state u and the parameter λ if we have observation data for ϕ at the end of the time interval t = T . We will minimize the cost function

J(u, λ) = inf w,v∈R J(w, v), (42)
where

J(u, λ) = α 2 u -u b 2 + 1 2
ϕ| t=T -ϕ obs 2 , α > 0, and ϕ is the solution to [START_REF] Zakharova | The new method of ARGO buoys system observation data interpolation[END_REF].

Thus, here we have

V 1 = α, V 2 = 0, V 3 = 1, Cϕ = ϕ| t=T .
In this case F ′ ϕ (ϕ, λ) = a, F ′ λ (ϕ, λ) = g, and the optimality system ( 11)-( 14) has the form:

   dϕ dt + aϕ = λg, t ∈ (0, T ) ϕ t=0 = u, (43)
   dϕ * dt -aϕ * = 0, t ∈ (0, T ) ϕ * t=T = ϕ obs -ϕ t=T , (44)
α(u -u b ) -ϕ * t=0 = 0, (45)
(g, ϕ * ) = T 0 g(t)ϕ * (t)dt = 0. (46)
It is easy to see that the problem of data assimilation ( 41)-(42) has a unique solution

λ = λ opt = ϕ obs -ϕ 0 u b ϕ 1 , u = u opt = u b , (47)
where

ϕ 0 = e -aT , ϕ 1 = T 0 e -a(T -t ′ ) g(t ′ )dt ′ .
Indeed, if u, λ have the form (47), the solution of problem [START_REF] Zakharova | The new method of ARGO buoys system observation data interpolation[END_REF] satisfies ϕ| t=T = ϕ obs , and the functional J from (42) attains its minimal value J = 0. In this case ϕ * = 0, and the optimality system (43)-( 46) is satisfied. Also, we will see below that the Hessian of J is positive definite, and it means the uniqueness of the solution u, λ.

Let us consider the response function in the form

G(ϕ, λ, u) = T 0 ϕ(t)dt. ( 48
)
Let a = 0. After assimilation, taking into account the solution of problem [START_REF] Zakharova | The new method of ARGO buoys system observation data interpolation[END_REF], we have

G(ϕ, λ, u) = u opt a (1 -e -aT ) + λ opt a   T 0 g(t)dt -ϕ 1   , (49)
where u opt , λ opt are given by (47). Then, by direct differentiation of G with respect to ϕ obs we have the gradient

dG dϕ obs = 1 aϕ 1   T 0 g(t)dt -ϕ 1   . ( 50
)
Let us now apply the algorithm ( 38)- [START_REF] Karagali | SST diurnal variability in the North Sea and the Baltic Sea, Remote Sens[END_REF] to compute the gradient of the function G. Since ∂G ∂ϕ = 1, (F ′ ϕ (ϕ, λ)) * = -a, then on the first step of the algorithm, we solve the problem (38) and get the solution

φ * (t) = 1 a (1 -e -a(T -t) ). (51)
Taking into account that ∂G/∂λ = ∂G/∂u = 0 and (F ′ λ (ϕ, λ)) * φ * = (g, φ * ), we get F = φ * (0), (g, φ * )

T , i.e., F = ( f0 , f ) T , where

f0 = φ * t=0 , f = T 0 g φ * dt = 1 a   T 0 g(t)dt -ϕ 1   . (52)
On the second step of the algorithm, one need to solve the equation HU = F with the Hessian H defined by the formulas (30)- [START_REF] Chavent | Local stability of the output least square parameter estimation technique[END_REF]. Since all the second order derivatives of F (ϕ, λ) equal zero, then it is easily seen that H in this case is defined by

HU = αw -φ * t=0 , - T 0 g(t)φ * (t)dt T , U = (w, v) T ,
where φ * is the solution of the adjoint problem

   dφ * dt -aφ * = 0, t ∈ (0, T ) φ * t=T = -φ t=T , (53)
and φ is the solution of the forward problem

   dφ dt + aφ = vg, t ∈ (0, T ) φ t=0 = w. (54)
Since φ t=T = wϕ 0 + vϕ 1 and

T 0 g(t)φ * (t)dt = -φ t=T T 0 e -a(T -t) g(t)dt = -φ t=T ϕ 1 , φ * t=0 = -φ t=T e -aT = -φ t=T ϕ 0 , we get HU = αw + ϕ 2 0 w + ϕ 0 ϕ 1 v, ϕ 0 ϕ 1 w + ϕ 2 1 v T , hence H is the 2 × 2 matrix H = α + ϕ 2 0 ϕ 0 ϕ 1 ϕ 0 ϕ 1 ϕ 2 1 . (55)
For α > 0, ϕ 0 , ϕ 1 = 0 the matrix H is positive definite, which confirms the existence and uniqueness of the solution to problem (42).

The solution of the system HU = F has the explicit form

w = - ϕ 0 f αϕ 1 + f0 α , v = (α + ϕ 2 0 ) αϕ 2 1 f - ϕ 0 αϕ 1 f0 . ( 56
)
On the third step of the algorithm, we need to solve problem [START_REF] Zalesny | Numerical modeling of marine circulation, pollution assessment and optimal ship routes[END_REF]. Since F ′ λ (ϕ, λ) = g, the solution of this problem for t = T has the form

P 2 t=T = wϕ 0 + vϕ 1 = - ϕ 2 0 f αϕ 1 + (α + ϕ 2 0 ) f αϕ 1 = f ϕ 1 .
Finally, using [START_REF] Karagali | SST diurnal variability in the North Sea and the Baltic Sea, Remote Sens[END_REF], we get the gradient of G with respect to ϕ obs :

dG dϕ obs = CP 2 = P 2 | t=T = f ϕ 1 . (57)
From ( 57) and (52) we have

dG dϕ obs = 1 aϕ 1   T 0 g(t)dt -ϕ 1   . (58)
Thus, the gradient obtained by the algorithm ( 38)-( 40) exactly coincides with the value of the gradient obtained in (50) by direct differentiation, which is the expected result.

For a numerical example, we consider the problem ( 41)-(42) and the response function G in the form (48) for a = 1, g(t) = 1, α = 10 -5 . The exact value of the gradient dG dϕ obs is defined by the formula (50). Easy to see that for a = 1, g(t) = 1 it has the explicit form:

dG dϕ obs = T 1 -e -T -1,
and it does not require a numerical integration. Approximate values of the gradient were obtained numerically with the help of the algorithm (38)- [START_REF] Karagali | SST diurnal variability in the North Sea and the Baltic Sea, Remote Sens[END_REF], where the problems [START_REF] Marchuk | Adjoint Equations and Perturbation Algorithms in Nonlinear Problems[END_REF], [START_REF] Zalesny | Numerical modeling of marine circulation, pollution assessment and optimal ship routes[END_REF] were solved using the simplest explicit scheme in time. Table 1 presents for comparison the exact and approximate values of the gradient for different lengths T of the assimilation window and for different time steps τ . One can see from the table that for each case the gradient values obtained by the considered algorithm coincide with the exact values of the gradient with the accuracy O(τ ). Besides, the gradient dG dϕ obs rises with the increase of T , therefore, the sensitivity of the response function with respect to observation errors is increasing, which is natural for a larger assimilation window T . Table 2. The experiment with different parameters α, T = 10, τ = 0.01 Note that in this example, the exact gradient defined by the formula (50) does not depend on α, however, the algorithm (38)-( 40) involve α as a parameter. Table 2 shows that the resulting approximate values of the gradient obtained by the algorithm also do not change with α.

Assimilation window

T = 1 T = 2 T = 5 T = 10 T =

Application: data assimilation problem for a sea thermodynamics model

We consider the sea thermodynamics problem in the form [START_REF] Marchuk | Mathematical Models in Geophysical Hydrodynamics and Numerical Methods for their Realization[END_REF]:

T t + ( Ū , Grad)T -Div(â T • Grad T ) = f T in D × (t 0 , t 1 ), T = T 0 for t = t 0 in D, -ν T ∂T ∂z = Q on Γ S × (t 0 , t 1 ), ∂T ∂n = 0 on Γ w,c × (t 0 , t 1 ), Ū (-) n T + ∂T ∂n = Q T on Γ w,op × (t 0 , t 1 ), ∂T ∂n = 0 on Γ H × (t 0 , t 1 ), (59)
where T = T (x, y, z, t) is an unknown temperature function, t ∈ (t 0 , t 1 ), (x, y, z) The other notations and a detailed description of the problem statement can be found in [START_REF] Agoshkov | Numerical algorithm of variational assimilation of the ocean surface temperature data[END_REF].

∈ D = Ω × (0, H), Ω ⊂ R 2 , H = H(x, y) is the function of the bottom releif, Q = Q(x, y, t) is the total heat flux, Ū = (u, v, w), a T = diag((a T ) ii ), (a T ) 11 = (a T ) 22 = µ T , (a T ) 33 = ν T , f T = f T (x,
Problem (59) can be written in the form of an operator equation:

T t + LT = F + BQ, t ∈ (t 0 , t 1 ), T = T 0 , t = t 0 , ( 60
)
where the equality is understood in the weak sense, namely,

(T t , T ) + (LT, T ) = F ( T ) + (BQ, T ) ∀ T ∈ W 1 2 (D), (61)
in this case L, F , B are defined by the following relations:

(LT, T ) ≡ D (-T Div( Ū T ))dD + Γw,op Ū (+) n T T dΓ + D a T Grad(T ) • Grad( T )dD, F ( T ) = Γw,op Q T T dΓ + D f T T dD, (T t , T ) = D T t T dD, (BQ, T ) = Ω Q T z=0 dΩ,
and the functions a T , Q T , f T , Q are such that equality (61) makes sense. The properties of the operator L were studied by [START_REF] Agoshkov | Numerical algorithm of variational assimilation of the ocean surface temperature data[END_REF].

Due to (61), the equation ( 60) is considered in Y = L 2 (t 0 , t 1 ; (W 1 2 (D)) * ), and the operator 59) is linear in T, Q, however, written in the form (60), it is a particular case of the original problem [START_REF] Lions | Contrôle Optimal des Systèmes Gouvernés par des Équations aux Dérivées Partielles[END_REF], and all the reasoning and the methodology presented in Sections 2-3 are easily transferred to the case of problem (60), understood in a weak sense.

B : L 2 (Ω × (t 0 , t 1 )) → Y maps the function Q ∈ L 2 (Ω × (t 0 , t 1 )) into the function BQ ∈ Y such that (BQ, T ) = Ω Q T z=0 dΩ, ∀ T ∈ W 1 2 (D). Problem (
We consider the data assimilation problem for the sea surface temperature (see [START_REF] Agoshkov | Numerical algorithm of variational assimilation of the ocean surface temperature data[END_REF]). Suppose that the functions Q ∈ L 2 (Ω × (t 0 , t 1 )) and T 0 ∈ L 2 (D) are unknown in problem (59). Let also T obs (x, y, t) ∈ L 2 (Ω × (t 0 , t 1 )) be the function on Ω obtained for t ∈ (t 0 , t 1 ) by processing the observation data, and this function in its physical sense is an approximation to the surface temperature function on Ω, i.e. to T z=0 . We admit the case when T obs is defined only on some subset of Ω × (t 0 , t 1 ) and denote the indicator (characteristic) function of this set by m 0 . For definiteness sake, we assume that T obs is zero outside this subset.

Consider the data assimilation problem for the surface temperature in the following form: find T 0 and Q such that

       T t + LT = F + BQ in D × (t 0 , t 1 ), T = T 0 , t = t 0 J(T 0 , Q) = inf w,v J(w, v), (62)
where

J(T 0 , Q) = α 2 t1 t0 Ω |Q-Q (0) | 2 dΩdt+ β 2 D |T 0 -T (0) | 2 dD + 1 2 t1 t0 Ω m 0 |T z=0 -T obs | 2 dΩdt, (63) and Q
(0) = Q (0) (x, y, t), T (0) = T (0) (x, y, z) are given functions, α, β = const > 0.
Lemma 2. For α, β > 0 the variational data assimilation problem (62) has a unique solution.

Proof. Let T n 0 , Q n be a sequence minimizing J(T 0 , Q), i.e. J(T n 0 , Q n ) → inf T0,Q J(T 0 , Q), n → ∞. Since J(T 0 , Q) ≥ α 2 t1 t0 Ω |Q-Q (0) | 2 dΩdt+ β 2 D |T 0 -T (0) | 2 dD, ∀T 0 ∈ L 2 (D), Q ∈ L 2 (Ω×(t 0 , t 1 )),
then, for α, β > 0, the sequence T n 0 , Q n is bounded:

T n 0 L2(D) ≤ const, Q n L2(Ω×(t0,t1)) ≤ const.
Hence, there exists a weakly convergent subsequence (we denote it also by T n 0 , Q n ). The Hilbert spaces L 2 (D) and L 2 (Ω × (t 0 , t 1 )) are weakly closed, therefore, there exist elements

T 0 ∈ L 2 (D), Q ∈ L 2 (Ω × (t 0 , t 1 )) such that T n 0 → T 0 weakly in L 2 (D), and Q n → Q weakly in L 2 (Ω × (t 0 , t 1 )), i.e. (T n 0 , p) L2(D) → (T 0 , p) L2(D) , (Q n , q) L2(Ω×(t0,t1)) → (Q, q) L2(Ω×(t0,t1)) ∀p ∈ L 2 (D), q ∈ L 2 (Ω × (t 0 , t 1 )
). Let T n and T be the solutions of problem (60) for T n 0 , Q n and T 0 , Q, respectively. Then, for the difference we have

(T n -T ) t + L(T n -T ) = B(Q n -Q), t ∈ (t 0 , t 1 ), T n -T = T n 0 -T 0 , t = t 0 . (64)
The solution to problem (64) continuously depends on the initial value T n 0 -T 0 and the flux Q n -Q (a priori estimates are valid in the corresponding functional spaces) [START_REF] Lions | Contrôle Optimal des Systèmes Gouvernés par des Équations aux Dérivées Partielles[END_REF], therefore, T n → T weakly in Y , and

T n | z=0 → T | z=0 weakly in L 2 (Ω × (t 0 , t 1 )). The functional S(•) =
• 2 is known [START_REF] Lions | Contrôle Optimal des Systèmes Gouvernés par des Équations aux Dérivées Partielles[END_REF] to be lower semi-continuous in the weak topology, then lim inf J(T n 0 , Q n ) ≥ J(T 0 , Q), and, therefore, inf

w,v J(w, v) ≥ J(T 0 , Q). Hence, inf w,v J(w, v) = J(T 0 , Q), that is, T 0 , Q
gets the minimum to the functional J. This proves the lemma.

The optimality system determining the solution of the formulated variational data assimilation problem according to the necessary condition gradJ = 0 has the form:

T t + LT = F + BQ in D × (t 0 , t 1 ), T = T 0 , t = t 0 , (65)
-(T * ) t + L * T * = Bm 0 (T obs -T ) in D × (t 0 , t 1 ), T * = 0, t = t 1 , (66)
α(Q -Q (0) ) -T * = 0 on Ω × (t 0 , t 1 ), (67)
β(T 0 -T (0) ) -T * t=0 = 0 in D, (68)
where L * is the operator adjoint to L.

Here the boundary-value function Q plays the role of λ from Section 2, ϕ = T , the operator F has the form F (T, Q) = -LT + BQ, and

F ′ T = -L, F ′ Q = B.
Since the operator F (T, Q) is linear in this case and F ′′ T T = F ′′ QQ = F ′′ QT = 0, the Hessian H acting on some U = (w, ψ) T , w ∈ L 2 (D), ψ ∈ L 2 (Ω×(t 0 , t 1 )) is defined by the successive solution of the following problems:

∂φ ∂t + Lφ = Bψ, t ∈ (t 0 , t 1 ) φ t=t0 = w, (69)
- ∂φ * ∂t + L * φ * = -Bm 0 φ, t ∈ (t 0 , t 1 ) φ * t=t1 = 0, (70)
HU = (βw -φ * t=0 , αψ -B * φ * ) T . (71)
To illustrate the above-presented theory, we consider the problem of sensitivity of functionals of the optimal solution Q to the observations T obs . Let us introduce the following response function:

G(T ) = t1 t0 dt Ω k(x, y, t)T (x, y, 0, t)dΩ, (72)
where k(x, y, t) is a weight function related to the temperature field on the sea surface z = 0.

For example, if we are interested in the mean temperature of a specific region of the sea ω for z = 0 in the interval t -τ ≤ t ≤ t, then as k we take the function

k(x, y, t) = 1 (τ mes ω) if (x, y) ∈ ω, t -τ ≤ t ≤ t 0 else, (73)
where mes ω denotes the area of the region ω. Thus, the functional (72) is written in the form:

G(T ) = 1 τ t t-τ dt 1 mes ω ω T (x, y, 0, t)dΩ . (74)
Formula (74) represents the mean temperature averaged over the time interval t -τ ≤ t ≤ t for a given region ω. The response functions of this type are of most interest in the theory of climate change ( [START_REF] Marchuk | Adjoint Equations and Analysis of Complex Systems[END_REF], [START_REF] Marchuk | Adjoint Equations and Perturbation Algorithms in Nonlinear Problems[END_REF]).

In our notations the functional (72) may be written as

G(T ) = t1 t0 (Bk, T )dt = (Bk, T ) Y , Y = L 2 (D × (t 0 , t 1 )).
We are interested in the sensitivity of the response function G(T ), obtained for T after data assimilation, with respect to the observation function T obs .

By definition, the sensitivity is given by the gradient of G with respect to T obs :

dG dT obs = ∂G ∂T ∂T ∂T obs . ( 75
)
Since ∂G ∂T = Bk, then according to the theory presented in Section 5, to compute the gradient (75) we need to perform the following steps:

1) For k defined by ( 73

The last formula allows us to estimate the sensitivity of the functionals related to the mean temperature after data assimilation, with respect to the observations on the sea surface.

For numerical experiments have used the three-dimensional numerical model of the Baltic Sea hydrothermodynamics developed at the INM RAS on the base of the splitting method [START_REF] Zalesny | Numerical modeling of marine circulation, pollution assessment and optimal ship routes[END_REF] and supplied with the assimilation procedure [START_REF] Agoshkov | Numerical algorithm of variational assimilation of the ocean surface temperature data[END_REF] for the surface temperature T obs with the aim to reconstruct the heat fluxes Q and the initial state T 0 .

The parameters of the considered domain of the Baltic Sea and its geographic coordinates can be described as follows: σ-grid is 336×394×25 (the latitude, longitude, and depth, respectively). The first point of the "grid C" [START_REF] Zalesny | Numerical modeling of marine circulation, pollution assessment and optimal ship routes[END_REF] has the coordinates 9.406 • E and 53.64 • N. The mesh sizes in x and y are constant and equal to 0.0625 and 0.03125 degrees. The time step is ∆t = 5 minutes. The assimilation procedure worked only during some time windows. To start the assimilation procedure, the function T (0) was taken as a model forecast for the previous time interval.

The Baltic Sea daily-averaged nighttime surface temperature data were used for T obs . These are the data of the Danish Meteorological Institute based on measurements of radiometers (AVHRR, AATSR and AMSRE) and spectroradiometers (SEVIRI and MODIS) [START_REF] Karagali | SST diurnal variability in the North Sea and the Baltic Sea, Remote Sens[END_REF]. Data interpolation algorithms were used [START_REF] Zakharova | The new method of ARGO buoys system observation data interpolation[END_REF] to convert observations on computational grid of the numerical model of the Baltic Sea thermodynamics. The mean climatic flux obtained from the NCEP (National Center for Environmental Prediction) reanalysis was taken for Q (0) .

Using the hydrothermodynamics model mentioned above, which is supplied with the assimilation procedure for the surface temperature T obs , we have performed calculations for the Baltic Sea area where the assimilation algorithm worked only at certain time moments t 0 ; in this case t 1 = t 0 + ∆t. The aim of the experiment was the numerical study of the sensitivity of functionals of the optimal solution Q, T 0 to observation errors in the interval (t 0 , t 1 ).

We use the discretize-then-optimize approach, and for numerical experiments all the presented equations are understood in a discrete form, as finite-dimensional analogues of the corresponding problems, obtained after approximation. This allows us to consider model equations as a perfect model, with no approximation errors.

Let us present some results of numerical experiments.

The calculation results for t 0 = 41 hours 40 minutes (500 time steps for the model) are presented in Fig. 1 showing the gradient of the response function G(T ) defined by (74) and related to the mean temperature after data assimilation, with respect to the observations on the sea surface, according to (76)-(78). Here ω = Ω, τ = ∆t, t = t 1 , α = β = 10 -5 . We can see the sub-areas (in red) in which the response function G(T ) is most sensitive to errors in the observations during assimilation. The largest values of the gradient of G(T ) correspond to the points x, y lying near the boundary of the domain. This result is confirmed by the direct computation of the response function G(T ) according to (74) obtained after assimilation, by introducing perturbations into the observation data T obs .

The above studies allow to determine the sea sub-areas in which the response function related to the optimal solution is most sensitive to errors in the observations during variational data assimilation.

Conclusions

Numerical algorithms are considered to study the sensitivity of functionals of the optimal solution of variational data assimilation problem aimed at the reconstruction of unknown parameters and initial state of the model. The optimal solution obtained as a result of assimilation depends on the observations that may contain uncertainties. Computing the gradient of the functionals with respect to observations reduces to the solution of a non-standard problem which is a coupled system involving direct and adjoint equations with mutually dependent variables. Solvability of the non-standard problem is related to the properties of the Hessian of the original cost function. An algorithm to compute the gradient of the response function is developed. Numerical example for variational data assimilation problem related to sea surface temperature for the Baltic Sea thermodynamics model demonstrates the result of the gradient computation of the response function associated with the mean surface temperature. The presented algorithm may be used to determine the sea sub-areas in which the response functions of the optimal solution are most sensitive to errors in the observations during variational data assimilation.

Appendix: Proof of Theorem 1

Consider the system of perturbations ( 16)- [START_REF] Alifanov | Extreme Methods for Solving Ill-posed Problems with Applications to Inverse Heat Transfer Problems[END_REF]. We have here 4 problems and 4 unknowns δu, δλ, δϕ, δϕ * . The function δϕ obs is included in the right-hand side of [START_REF] Shutyaev | Sensitivity with respect to observations in variational data assimilation[END_REF], therefore, all the unknowns δu, δλ, δϕ, δϕ * will depend on δϕ obs . The expression [START_REF] Sun | Inverse Problems in Groundwater Modeling[END_REF] involves δϕ obs in the left-hand side, and we would like to represent the right-hand side of (20) through δϕ obs also, to obtain the explicit formula for the gradient dG dϕ obs . Let us introduce four adjoint variables P 1 ∈ Y , P 2 ∈ Y , P 3 ∈ Y p and P 4 ∈ X. By taking the inner product of ( 16) by P 1 , (17) by P 2 , (19) by P 3 and of ( 18) by P 4 and adding them, we get:

We would like the first three summands in the left-hand side of (80) be equal to the rihgt-hand side of [START_REF] Sun | Inverse Problems in Groundwater Modeling[END_REF], keeping the others summands to be zero, therefore, we put , and due to [START_REF] Sun | Inverse Problems in Groundwater Modeling[END_REF] the gradient of G is given by [START_REF] Schirber | Parameter estimation using data assimilation in an atmospheric general circulation model: From a perfect toward the real world[END_REF]. The theorem is proved.

-

  y, z, t) are given functions. The boundary of the domain Γ ≡ ∂D is represented as a union of four disjoint parts Γ S , Γ w,op , Γ w,c , Γ H , where Γ S = Ω (the unperturbed sea surface), Γ w,op is the liquid (open) part of vertical lateral boundary, Γ w,c is the solid part of the vertical lateral boundary, Γ H is the sea bottom, Ū (-) n = (| Ūn |-Ūn )/2, and Ūn is the normal component of Ū .

  ∂t + L * φ * = Bk, t ∈ (t 0 , t 1 ) put Φ = ( φ * t=0 , B * φ * ) T .2) Find U = (w, v) T by solving HU = Φ with the Hessian defined by (69)-(71).3) Solve the direct problem∂P 2 ∂t + LP 2 = Bv, t ∈ (t 0 , t 1 ) P 2 t=t0 = w.(77)4) Compute the gradient of the response function as dG dT obs = m 0 P 2 z=0 .

Figure 1 :

 1 Figure 1: The gradient of the response function G(T )

2 Y+ V 2 3 Yp+ V 1 4 X = 0 . 2 Y++

 2231402 ϕ (ϕ, λ)δϕ-F ′ λ (ϕ, λ)δλ, P 1 Y + -∂δϕ * ∂t -(F ′ ϕ (ϕ, λ)) * δϕ * -(F ′′ ϕϕ (ϕ, λ)δϕ) * ϕ * --(F ′′ ϕλ (ϕ, λ)δλ) * ϕ * + C * V 3 (Cδϕ -δϕ obs ), P δλ -(F ′′ λϕ (ϕ, λ)δϕ) * ϕ * --(F ′′ λλ (ϕ, λ)δλ) * ϕ * -(F ′ λ (ϕ,λ)) * δϕ * , P δu -δϕ * t=0 , P Using integration by parts and adjoint operators, we obtain δϕ, -∂P 1 ∂t -(F ′ ϕ (ϕ, λ)) * P 1 -(F ′′ ϕϕ (ϕ, λ)P 2 ) * ϕ * -(F ′′ λϕ (ϕ, λ)P 3 ) * ϕ * + C * V 3 CP δϕ t=T , P

Table 1 .

 1 100 The experiment with different assimilation windows, α = 10 -5

	Exact gradient	0.5819 1.313 4.034	9	99
	Approximate
	gradient, τ = 0.1 0.5193 1.246 3.978	8.95	98.949
	Approximate
	gradient, τ = 0.01 0.5815 1.302 4.027	8.995	98.997
	Parameter α	α = 10 -5 α = 10 -2 α = 1
	Exact gradient	9	9	9
	Approximate
	gradient	8.995	8.995	8.995

  Yp-δϕ obs , V 3 CP 2 ′′ ϕϕ (ϕ, λ)P 2 ) * ϕ * -(F ′′ λϕ (ϕ, λ)P 3 ) * ϕ * + C * V 3 CP 2 , δϕ V 2 P 3 -(F ′′ ϕλ (ϕ, λ)P 2 ) * ϕ * -(F ′′ λλ (ϕ, λ)P 3 ) * ϕ * -(F ′ λ (ϕ, λ)) * P 1 , δλ

				+ P 2 t=0 -P 4 , δϕ *	t=0	X	= V 3 CP 2 , δϕ obs	Y obs	.
				1 t=T	X	-δu, P 1 t=0	X	+ δϕ * ,	∂P 2 ∂t	-F ′ ϕ (ϕ, λ)P 2 -F ′ λ (ϕ, λ)P 3	Y	+
	+ δϕ *	t=0 , P 2 t=0
						Y obs	+ δu, V 1 P 4	X	-δϕ *	t=0 , P 4	X	= 0.	(79)
	Hence,
	-	∂P 1 ∂t	-(F ′ ϕ (ϕ, λ)) Yp	+
	+ V 1 P 4 -P 1 | t=0 , δu	X	+ P 1 t=T , δϕ t=T		X	+	∂P 2 ∂t	-F ′ ϕ (ϕ, λ)P 2 -F ′ λ (ϕ, λ)P 3 , δϕ *	Y	+

X + δλ, V 2 P 3 -(F ′′ ϕλ (ϕ, λ)P 2 ) * ϕ * -(F ′′ λλ (ϕ, λ)P 3 ) * ϕ * -(F ′ λ (ϕ, λ)) * P 1 * P 1 -(F Y + +

  ∂P 1 ∂t -(F ′ ϕ (ϕ, λ)) * P 1 -(F ′′ ϕϕ (ϕ, λ)P 2 ) * ϕ * -(F ′′ λϕ (ϕ, λ)P 3 ) * ϕ * + C * V 3 CP 2 = P 3 -(F ′′ ϕλ (ϕ, λ)P 2 ) * ϕ * -(F ′′ λλ (ϕ, λ)P 3 ) * ϕ * -(F ′ λ (ϕ, λ)) * P 1 = ′ ϕ (ϕ, λ)P 2 -F ′ λ (ϕ, λ)P 3 = 0, P 2 t=0 -P 4 = 0.Thus, if P 1 , P 2 , P 3 , P 4 are the solutions of (21)-(24), we get from from (80):

								∂G ∂ϕ	,
			V 1 P 4 -P 1 t=0 =	∂G ∂u	,
	and
	V 2 ∂G ∂λ	, P 1 t=T = 0,
	∂P 2 ∂t -F ∂G ∂ϕ , δϕ Y +	∂G ∂λ	, δλ	Yp	+	∂G ∂u	, δu

X = V 3 CP 2 , δϕ obs Y obs

Acknowledgments

This work was carried out within Russian Science Foundation project 19-71-20035 (studies in Sections 2-3), and the project 18-01-00267 of the Russian Foundation for the Basic Research. The authors are greatly thankful to the reviewers for providing valuable comments that helped to improve the paper.