Minos Axenides
email: axenides@inp.demokritos.gr

Emmanuel Floratos
email: mflorato@phys.uoa.gr

Stam Nicolis
email: stam.nicolis@lmpt.univ-tours.fr

Exponential mixing of all orders for Arnol'd cat map lattices

We show that the recently introduced classical Arnol'd cat map lattice field theories, which are chaotic, are exponentially mixing to all orders. Their mixing times are welldefined and are expressed in terms of the Lyapunov exponents, more precisely by the combination that defines the inverse of the Kolmogorov-Sinai entropy of these systems. We prove by an explicit recursive construction of correlation functions, that these exhibit l-fold mixing for any l = 3, 4, 5, . . .. This computation is relevant for Rokhlin's conjecture, which states that 2-fold mixing induces l-fold mixing for any l > 2. Our results show that 2-fold exponential mixing, while being necessary for any l-fold mixing to hold it is nevertheless not sufficient for Arnol'd cat map lattice field theories. The correspondence principle implies that these mixing times, also, control the scrambling of the underlying quantum system for short times.

Introduction and Motivation

Fast information scrambling has occupied center stage in theoretical physics research in recent times. It forms an integral part of the so-called information loss paradox in quantum black holes and it is believed to be relevant for its long sought-after unitary resolution, thereby resolving the incompleteness of general relativity, using the principles of quantum mechanics [START_REF] Bousso | Snowmass White Paper: Quantum Aspects of Black Holes and the Emergence of Spacetime[END_REF].

Controversial issues of locality and unitarity between quantum theory and black hole physics have given rise to new concepts and principles in quantum gravity such as those of the holographic and complementarity principles and AdS/CFT correspondence, with black holes having been conjectured to be the fastest information scramblers in nature [START_REF] Hayden | Black holes as mirrors: Quantum information in random subsystems[END_REF][START_REF] Sekino | Fast Scramblers[END_REF]. In essence they highlight that the no-cloning theorem of quantum information imposes constraints on the speed of information spreading within infalling matter on the black hole, due to the interaction with its horizon degrees of freedom. Its characteristic timescale has been estimated to be proportional to the logarithm of the entropy of the black hole. It is called the scrambling time, a conjectured lower bound to the time of information propagation in the universe, with the search for a microscopic many-body system that can saturate it being an area of topical research.

More specifically the ingredients of regular local field theories have been shown to be inadequate to accommodate fast information processing among the near horizon black hole microscopic degrees of freedom; this has led to the study of how nonlocality can be taken into account [START_REF] Hayden | Black holes as mirrors: Quantum information in random subsystems[END_REF][START_REF] Almheiri | Black Holes: Complementarity or Firewalls?[END_REF][START_REF] Harlow | TF1 Snowmass Report: Quantum gravity, string theory, and black holes[END_REF]. A further property, that has emerged as relevant is that of chaoticity. This has, in turn, motivated the study of chaos in matrix models [START_REF] Iizuka | Matrix Models for the Black Hole Information Paradox[END_REF], conformal field theories and a large class of SYK type of models [START_REF] Maldacena | Remarks on the Sachdev-Ye-Kitaev model[END_REF] which have been dubbed chaotic field theories [START_REF] Cvitanovic | Chaotic field theory: A Sketch[END_REF].

In a recent work we have studied the AdS 2 near horizon geometry of extremal black holes by means of a modular discretization thereof [START_REF] Axenides | Modular discretization of the AdS 2 /CFT 1 holography[END_REF]. This discretization is relevant for probing the non-local features of the geometry using single particle excitations. Among such excitations the Arnol'd cat map stands out, since its isometries are those of the near horizon geometry and its dynamics can capture fast scrambling in a way that saturates the fast scrambling bound [START_REF] Axenides | Chaotic Information Processing by Extremal Black Holes[END_REF][START_REF] Axenides | The quantum cat map on the modular discretization of extremal black hole horizons[END_REF]. Furthermore, the model has been shown to admit a well defined continuum limit [START_REF] Axenides | The arithmetic geometry of AdS 2 and its continuum limit[END_REF] showing how a discrete short distance near horizon spacetime can underlie an emergent continuum long distance classical limit.

In followup work we have taken a first step in implementing the conditions of nonlocality and strong chaos for many-body systems through the construction of lattice field theories, in which we replace the standard harmonic oscillator constituents with chaotic interacting constituents, namely, the so called Arnol'd cat map [START_REF] Axenides | Arnol'd cat map lattices[END_REF][START_REF] Gutkin | Linear encoding of the spatiotemporal cat[END_REF]. The resulting evolution matrix of the n-body system is represented by elements of the symplectic group Sp 2n [Z] . The model has been shown to exhibit strong classical chaos as well as non-locality, which leaves a direct imprint on the (spatiotemporal) chaotic benchmark of the model, its Kolmogorov-Sinai entropy.

In the present work we pursue the study of the chaotic properties of the many-body system, by elucidating its mixing properties, the classical avatars of scrambling. We will show that it exhibits fast exponential mixing of all orders in the number of observables. To this end we adapt an operator method developed by de Bièvre [START_REF] Bievre | Chaos, quantization and the classical limit on the torus[END_REF], firstly for the case of two observable functions in the computation of the decay rates of their correlation function. Moreover we generalize his method for any number of observables for the Arnol'd cat map lattice, deducing that it is a strongly l-fold mixing system, for any l = 2, 3, . . . and that the mixng time is 1/S K-S ,, where S K-S is the Kolmogorov-Sinai entropy.

The role of deterministic chaos [START_REF] Zaslavsky | Chaos in dynamic systems[END_REF] in rendering mixing phenomena fast and "efficient" is well recognized and can be understood within the so-called ergodic hierarchy classification of hyperbolic dynamical systems [START_REF] Berkovitz | The ergodic hierarchy, randomness and hamiltonian chaos[END_REF]. Mixing is the classical analog of the fast scrambling of information, which is the rapid spreading of an initially localized perturbation, as it evolves towards a homogeneous stationary state.

Our results pertain, also, to the well known Rokhlin conjecture [START_REF] Rokhlin | On endomorphisms of compact commutative groups[END_REF][START_REF] Kosenko | On Rokhlin's multiple mixing problem[END_REF] of whether 2-fold mixing implies multiple l-fold mixing and more generally mixing of all higher orders. We show, by explicit calculation, that 2-fold mixing is necessary but not sufficient for any l-fold (l > 2) mixing to hold. We proceed, to this end, in steps, by establishing, first, both the 2fold and l-fold mixing property of a single ACM. Subsequently we on the case of the extended n-body coupled Arnold cat map systems through the explicit computation of l-fold correlation functions. We find, by explicit calculation, that the correlation functions show exponential decay behavior and identify the mixing time as the inverse of the Kolmogorov-Sinai entropy. The reason that 2-fold mixing isn't sufficient for ensuring l-fold mixing can be understood from the recursive construction of the l + 1st correlation function, which is the sum of two terms, the first being the l-th correlation function, while the second is a remainder term that can be shown, also, to tend to zero for large times-with subleading behavior.

The plan of the paper is the following:

In section 2 we review the salient features of our previous paper [START_REF] Axenides | Arnol'd cat map lattices[END_REF], namely the construction of the evolution operator for n maps, the analytical computation of the spectrum of the Lyapunov exponents, from which we have obtained the Kolmogorov-Sinai entropy, S K-S .

In sections 3 4 we discuss the definition of simple and multiple mixing of chaotic dynamical systems and we provide an explicit calculation of the corresponding mixing times for the case of a single Arnol'd cat map.

In section 5 we extend the previous calculation to the case of n, symplectically interacting, Arnol'd cat maps, for any n. We find that the corresponding, simple, mixing time is equal to 1/S K-S .

Our conclusions and discussion for directions of further inquiry are set forth in section 6.

2 Interacting Arnol'd cat maps from symplectic couplings of n, k-Fibonacci sequences

In this section we review the salient features of our previous paper [START_REF] Axenides | Arnol'd cat map lattices[END_REF], pertaining to the construction of the evolution operator for n interacting Arnol'd cat maps and the analytical computation of the spectrum of the Lyapunov exponents, from which we have obtained the Kolmogorov-Sinai entropy, S K-S .

We consider a dynamical system of n degrees of freedom, whose space of states is the torus T 2n , the 2n-dimensional torus of radius 1. Therefore any state can be identified with a point x ∈ T 2n .

We shall describe the time evolution of this system by successive applications on a point x ∈ T 2n of elements, M, of the symplectic group over the integers, Sp 2n [Z]. The evolution is, thus, discrete in time. In equations,

x m+1 = x m M mod 1 (2.1)
where x m = (q m , p m ) describes the state of the system at the m-th time tick. Here the positions and momenta take values in T n . These equations describe the evolution of n "particles", from the tick m to the tick m + 1.

Since M ∈ Sp 2n [Z], it preserves the symplectic structure,

J = 0 -I I 0 (2.2)
in the sense that

M T JM = J (2.3)
Upon decomposing M into four blocks of n × n matrices

M = A B C D (2.4)
we obtain the constraints on the blocks,

A T D -C T B = I n×n A T C = C T A B T D = D T B (2.5)
It is known that Sp 2n [Z] is generated by the elements

M R = I R 0 I , M L = I 0 L I and D S = S T 0 0 S -1 (2.6)
where R and L are symmetric matrices; all matrices have integer entries.

In what follows we are interested in those evolution matrices, M, that have strictly positive eigenvalues-the reason is that this property is sufficient for the dynamics to be chaotic. Such symplectic matrices are called hyperbolic and the corresponding dynamical systems are called hyperbolic or Anosov.

It can be shown that the symplectic property implies that the eigenvalues come in pairs, (ρ

(i) + , ρ (i) -= 1/ρ (i)
+ < 1). For hyperbolic matrices this implies the existence of n planes, spanned by the corresponding eigenvectors: The dynamics is expanding along the eigenvectors, v We shall focus here on evolution matrices M that are hyperbolic and symmetric, so the eigenvectors v (i) ± span orthogonal eigenspaces and v

(i) + and v (i) -are orthogonal for each i = 1, 2, . . . , n.
This class of matrices includes, in particular, the generalization of the so-called Arnol'd cat map, defined by the 2 × 2 matrix

M = 1 1 1 2 (2.7)
to the case of n such maps, in interaction. Moreover, since the phase space is compact, this class of maps exhibits strong mixing as well as l-fold mixing of any order, l, as we shall show in the following. To this end we shall use the Lyapunov exponents and the Kolmogorov-Sinai entropy, that were calculated in closed form in [START_REF] Axenides | Arnol'd cat map lattices[END_REF].

In order to study the effects of coupling of n Arnol'd cat maps, each of which is defined on a lattice of n sites, we associate to each site a two-dimensional torus, with dynamics described by a single Arnol'd cat map. The total phase space of the system will be T 2n , the 2n-dimensional torus and the proposed dynamics will be described by appropriate elements of the symplectic group, Sp 2n [Z]. The 2n-dimensional symplectic maps will allow couplings of various degrees of locality and strength.

To begin with, we shall show that these maps can be constructed as iteration matrices of n coupled Fibonacci sequences.

We start from the relation between the ACM and the Fibonacci sequence of integers.

The Fibonacci sequence is one of the integer sequences, which has been studied, for a long time and there are journals dedicated to its properties and their applications.

The definition is given by the relations

f 0 = 0; f 1 = 1 f m+1 = f m + f m-1 (2.8)
which can be written in matrix form

f m f m+1 = 0 1 1 1 A f m-1 f m (2.9)
The matrix A is not a symplectic matrix, but it satisfies

A T JA = -J (2.10) for n = 1.
We remark that the Arnol'd cat map, acting on the torus T 2 , can be written as

1 1 1 2 = A 2 (2.11)
Eq. (2.10) implies that [A 2 ] T JA 2 = J, therefore that A 2 is symplectic. It's possible to generalize the Fibonacci sequence in the following way:

g m+1 = kg m + g m-1 (2.12)
with g 0 = 0 and g 1 = 1 and k is a positive integer. This is known as the "k-Fibonacci" sequence [START_REF] Horadam | A Generalized Fibonacci Sequence[END_REF]. We may solve eq. (2.12) by g m ≡ Cρ m . The characteristic equation for ρ reads

ρ 2 -kρ -1 = 0 ⇔ ρ ± (k) = k ± √ k 2 + 4 2 (2.13)
and express g m as a linear combination of the ρ ± , upon taking into account the initial conditions:

g m = A + ρ + (k) m + A -ρ -(k) m = ρ + (k) m -(-) m ρ + (k) -m √ k 2 + 4 (2.14)
In matrix form

g m g m+1 = 0 1 1 k A(k) g m-1 g m (2.15)
Similarly as for the usual Fibonacci sequence, we may show, by induction, that

A(k) m = g m-1 g m g m g m+1 (2.16)
We remark that det A(k) m = (-) m = g m+1 g m-1 -g 2 m , for m = 1, 2, . . . and that lim m→∞ (g m+1 /g m ) = ρ + (k), which, for k = 1, is the golden ratio, for k = 2 is the silver ratio and, for k > 2 are generalizations thereof. 2 is the greatest Lyapunov exponent of the generalized evolution operator, i.e. the k-Arnol'd cat map, A(k) 2 and is an increasing function of k.

λ + (k) ≡ log ρ + (k)
Coupled Fibonacci sequences have been considered in the literature, for instance in [START_REF] Rathore | Generalized Coupled Fibonacci Sequences[END_REF]. However, in these papers the possible applications to Hamiltonian dynamics, were not the topic of interest and moreover the corresponding maps were not symplectic.

After this review of the single Arnol'd cat map, we proceed to the study of how many such maps can interact, in a way such that the evolution matrix is an element of Sp 2n [Z]. To this end, we shall, once more, use the correspondence between the Arnol'd cat maps and the Fibonacci sequences.

We start by coupling n = 2 Fibonacci sequences, {f m } and {g m } (but we write the expressions in a way that generalizes immediately to arbitrary n):

f m+1 = a 1 f m + b 1 f m-1 + c 1 g m + d 1 g m-1 g m+1 = a 2 g m + b 2 g m-1 + c 2 f m + d 2 f m-1 (2.17)
where the a i , b i , c i , d i , i = 1, 2 are integers, f 0 = 0 = g 0 and f 1 = 1 = g 1 are the initial conditions and m = 1, 2, 3, . . .. In matrix form, these read

X m+1 ≡     f m g m f m+1 g m+1     =     0 0 1 0 0 0 0 1 b 1 d 1 a 1 c 1 d 2 b 2 c 2 a 2         f m-1 g m-1 f m g m     Xm (2.18)
Let us define the 2×2 matrices

D ≡ b 1 d 1 d 2 b 2 C ≡ a 1 c 1 c 2 a 2 (2.19)
in terms of which the one-time-step evolution equation (2.18) can be written in block form as

X m+1 = 0 n×n I n×n D C X m (2.20)
In analogy with the case of a single Fibonacci sequence and its relation with the Arnol'd cat map, we impose the constraint (cf. (2.10))

0 n×n I n×n D C T J 0 n×n I n×n D C = -J (2.21)
This condition implies that

D = I n×n C = C T (2.22) Therefore a 1 = k 1 , a 2 = k 2 , c 1 = c 2 = c.
This implies, in particular, that the coupling between the sequences is the same for both, in order for the square of the evolution matrix to be symplectic. As we shall explain below, the generalization of this property is that, for more sequences, the corresponding coupling matrix must be symmetric. This is important in order to define chains of interacting k-Arnol'd cat maps.

In terms of these parameters, the recursion relations take the form

f m+1 = k 1 f m + f m-1 + cg m g m+1 = k 2 g m + g m-1 + cf m (2.23)
and can be identified as describing a particular coupling between a k 1 -and a k 2 -Fibonacci sequence. This particular coupling is determined by the condition that the square of the evolution matrix is an element of Sp 4 [Z]:

A = 0 1 1 C ⇒ M = A 2 = 1 C C 1 + C 2 (2.24)
The generalization to a chain of n k-Fibonacci sequences, with nearest-neighbor couplings, proceeds as follows: We choose two diagonal matrices, of positive integers,

K IJ = K I δ IJ and G IJ = G I δ IJ , with I, J = 1, 2, . . . , n.
We now define the translation operator, P along a closed chain, by P I,J = δ I-1,J mod n. The periodicity is expressed by the fact that P n = I n×n . Moreover, P is orthogonal, since PP T = I n×n . Now we can define the coupling matrix for n sequences as

C = K + PG + GP T (2.25)
The corresponding 2n × 2n evolution matrix, A is given by

A = 0 n×n I n×n I n×n C (2.26)
and satisfies the relation

A T JA = -J. Its square, M = A 2 = I n×n C C I n×n + C 2 (2.27) therefore satisfies the relation M T JM = J, showing that M ∈ Sp 2n [Z].
Since A is symmetric, (from the property that C = C T ), M is positive definite and its eigenvalues come in pairs, (ρ, 1/ρ), with ρ > 1 (and the corresponding eigenvectors are orthogonal). This property implies that, for all matrices K and G this system of coupled maps is hyperbolic. An important special case arises if we impose translation invariance along the chain, i.e. K I = K and G I = G for all I = 1, 2, . . . , n.

Let us now consider the case of the open chain. The only change involves the operator P, which, now, must be defined as P IJ = δ I-1,J , for I, J = 1, 2, . . . , n. Due to the absence of the mod n operation, the "far non-diagonal" (upper right and lower left) elements are, now, zero. This express the property that the n-th Fibonacci is not coupled to the first one (and vice versa).

For both, closed or open, chains, we observe certain algebraic properties of the evolution matrix, A.

The k-Fibonacci sequence has the important property that the elements of the matrix A(k) m are arranged in columns of consecutive pairs of the sequence. We shall show that this property can be generalized for n interacting k-Fibonacci sequences as follows:

Theorem 1. The m-th power of the evolution matrix, A (cf. eq. ( 2 Having constructed a large class of symplectic many-body maps, that describe the dynamics of n k-Arnol'd cat maps, we wish to understand their chaotic behavior, as they act on any initial condition of T 2n , for n > 1.

.26)) can be written as

A m = C m-1 C m C m C m+1 (
These particular symplectic matrices generalize the chaotic behavior of one Arnol'd cat map, acting on T 2 , to that of n maps, acting on T 2n . Indeed, starting with initial conditions x 0 ∈ T 2n = 0, which have irrational components, the evolution matrix M, as it acts on x 0 , defines an orbit,x m = x 0 M m , that will, in the limit m → ∞, cover the whole torus T 2n ; the full phase space is the attractor of the map. This means that the map is ergodic.

Furthermore, this map has the, additional, property of being strongly chaotic, which means that it has positive Kolmogorov-Sinai entropy. It is possible to tune the parameters K and G (in the translation invariant case) so that no Lyapunov exponent is equal to zero (which would imply the existence of a conservation law). In this cse the map is maximally hyperbolic, i.e. it is an "Anosov C-system".

Another quantitative measure of the "long time" chaotic behavior of the orbits of a hyperbolic map is provided by the properties of the time correlation functions of "observables", i.e. functions on the phase space, T 2n . These properties define the mixing behavior of the map and can be expressed by the fact that the connected correlation functions decay to zero, for long times. This, thus, raises the problem of computing them and, in particular whether they decay exponentially, thereby defining the spectrum of mixing times.

On the other hand, another way to look at the chaotic properties of a dynamical system, is to consider, if possible, the full set of its unstable periodic orbits. For large periods, the corresponding periodic orbit should approach (this is called "shadowing") the chaotic orbits, which fill the phase space [START_REF] Auerbach | Exploring chaotic motion through periodic orbits[END_REF]. Fortunately, for our system, of n interacting Arnol'd cat maps, the full set of unstable periodic orbits is produced by all the rational points of T 2n , taken as initial conditions. So the problem reduces to finding periodic orbits, with "very large" periods. This problem is difficult, because the periods are random functions of the-commondenominator, N, of the rational points, that are initial conditions [START_REF] Axenides | Arnol'd cat map lattices[END_REF].

In the following section we shall, thus, introduce the notion of mixing and the mixing time for ergodic systems and we shall present the method for obtaining a bound on the mixing time for the case of the single cat map; then we shall do the same calculation for two, symplectically coupled, cat maps, to see what are the issues that arise. Finally, we shall present the generalization to n symplectically coupled maps, where, in addition to the dependence on the coupling, the issue of the range becomes of interest.

The mixing time for one cat map

Mixing is the ubiquitous phenomenon of blending together distinct many body matter systems from an initial inhomogeneous state to a final homogeneous uniform configuration [START_REF] Cornfeld | Ergodic theory[END_REF]. The mixing time is the characteristic transient period it takes for an initial local perturbation to delocalize and spread in a many-body system attaining a uniform and homogenized final state. What sets the scale of mixing time as well as its precise determination is both a conceptual and a computational challenge.

The development of mathematical methods associated with the Ergodic Theory have been brought into prominence [START_REF] Cornfeld | Ergodic theory[END_REF][START_REF] Sinai | Introduction to ergodic theory[END_REF] in the study of "Ergodic Mixing" through the study of stochasticity in measure-preserving dynamical systems.

For a discrete-time dynamical system, T : M → M , which preserves a probability measure µ a "strong mixing"' condition can be formulated for two sets of points A and B on a constant energy surface E , as follows:

µ (A ∩ T n (B)) n-→∞ ----→ µ (A) µ (B) (3.1)
Intuitively the mixing condition states that a dynamical system is strongly mixing whenever any two observable events A and B which occur at separate time instances specified by the action of T n on B become independent in the infinite time separation limit n → ∞.

Equivalently its standard diagnostic quantifier is the decay of the correlation function for any pair of observable functions f, g :

M → C C n (f, g * ) = f • T n • g * dµ - f dµ • g * dµ n-→∞ ----→ 0 (3.2)
The correlation function C n (f, g * ) or the self-correlation for a single function C n (f, f * ) fall off to zero either polynomially, exponentially (rapid mixing) or even super exponentially in the long time limit n → ∞ [START_REF] Baladi | Decay of correlations[END_REF][START_REF] Pollicott | Exponential Mixing: Lectures from Mumbai[END_REF].

We shall show explicitly that the Arnol'd cat map on T 2 is exponentially mixing and that its mixing time is given by the expression

τ mix (Arnol ′ d) = 1 log ρ + (3.3)
where ρ + = (3 + √ 5)/2. This is, in fact equal to the inverse of the Kolmogorov-Sinai entropy of the system. This result is known and in the next section we shall generalize it to the case of symplectically coupled Arnol'd cat maps.

We start with some preparatory material: Any smooth and square integrable observable on T 2 has a uniformly convergent Fourier series

f (x) = k∈Z×Z c k e 2πik•x (3.4)
The set of all these observables defines a Hilbert space, H(T 2 ), with inner product

f, g = T 2 d 2 x f (x)g * (x) = k∈Z×Z c k • d * k (3.5) where g(x) = k∈Z×Z d k e 2πik•x (3.6)
and the norm is defined by

||f || 2 = T 2 d 2 x f (x)f * (x) = k∈Z×Z |c k | 2 (3.7)
For any f, g ∈ H(T 2 ) we have the Cauchy-Schwarz inequality

| f, g | ≤ ||f ||||g|| (3.8)
With these standard preliminaries, we proceed with the evaluation of the correlation functions for the Arnol'd cat map.

C n (f, g * ) = d 2 x f (T n x)g * (x) - d 2 x f (x) d 2 x g * (x) = f • T n , g * -f g * (3.9)
We follow the procedure sketched in ref. [START_REF] Bievre | Chaos, quantization and the classical limit on the torus[END_REF]. The idea is to choose as functions the eigenfunctions of a particular operator and show that, for these,

|C n (f, g * )| ≤ const × e -n/τ .
Let us recall the argument: The ACM has eigenvalues

ρ ± = 3 ± √ 5 2 (3.10)
and corresponding eigenvectors,

u ± = 1 √ 1 + ρ ± 1 ρ ± -1 (3.11)
Au ± = ρ ± u ± . Since ρ + ρ -= 1 and they're real (since the matrix is symmetric) ρ + > 1 and

ρ -= 1/ρ + < 1.
Now, let us define the operators

D ± ≡ - i 2π u ±,I δ IJ ∂ x J (3.12)
Since the torus, T 2 , is a compact manifold, the spectrum of the operators D ± is discrete and can be labeled by two integers, n I . Indeed, it is straightforward to check that the functions

e k (x) = e 2πik•x (3.13)
with k ∈ Z × Z and x ∈ T 2 are eigenfunctions of D ± :

D ± e k (x) = k • u ± e k (x) (3.14)
Since the components of the u ± are irrational numbers, k • u ± = 0 for any k ∈ Z × Z, therefore the inverses, [D ± ] -1 ,

[D ± ] -1 e k (x) = 1 k • u ± e k (x) (3.15)
are well-defined. We proceed below with the details of the evaluation of the correlation function C n (f, g * ) in a more explicit form.

First we split the Fourier sums of the functions f and g as follows:

C n (f, g * ) = k,l∈Z×Z c k d * l T 2 d 2 x e 2πik•T n x e -2πil•x -c 0 d * 0 = k =(0,0) c k d * 0 T d 2 x e 2πik•T n x + l =(0,0) c 0 d * l T 2 d 2 x e -2πil•x + k =(0,0),l =(0,0) c k d * l T 2 d 2 x e 2πik•T n x e -2πil•x
(3.16) The first two terms, that involve sums over k = (0, 0) or l = (0, 0) are zero (since the integrals are δ-functions on k = (0, 0) and l = (0, 0) respectively), therefore only the last term survives:

C n (f, g * ) = k =(0,0),l =(0,0) c k d * l T 2 d 2 x e 2πik•T n x e -2πil•x (3.17)
Therefore, the operators

D ± ≡ -(i/(2π))u ± • ∇ x act on the function f (x) = e 2πik•(A n x) as D ± e 2πik•(A n x) = u ± • kρ n ± e 2πik•(A n x)
as can be checked by direct calculation. These calculations now lead to the following statement:

Proposition 1. T 2 d 2 x D ± f (A n x) D -1 ± g * (x) = - T 2 d 2 x f (A n x) g * (x)
which can be proved using integration by parts.

Since any function on the torus can be expanded in plane waves, we deduce that the correlation function of any two functions on the torus, that are sufficiently smooth for their Fourier expansions (and those of their derivatives) to converge, will show mixing behavior with the same mixing time.

Using Proposition 1, this can be established as follows:

C n (f, g * ) = λ n ± T 2 d 2 x   k =(0,0) c k k • u ± e 2πik•T n x     l =(0,0) d * k 1 l • u ± e -2πil•x   (3.18)
We now apply the Cauchy-Schwarz inequality:

| f, g * | 2 ≤ ||f || 2 ||g * || 2 (3.19)
to deduce that

|C n (f, g * )| 2 ≤ ρ 2n ± T 2 d 2 x k =(0,0) c k k • u ± e 2πik•x 2 T 2 d 2 x l =(0,0) d * k l • u ± e -2πil•x 2 (3.20)
The integrals do not depend on n, so the best bound on how

|C n (f, g * )| vanishes, as n → ∞, is obtained in the form |C n (f, g * )| ∼ e -n log ρ + = e -nλ + (3.21)
as n → ∞.

Multiple mixing for the single Arnol'd cat map

We now proceed to extend the notion of mixing to any number of observables in both the ergodic theory of measure preserving transformations as well as to their corresponding diagnostic quantifiers of decay of higher order correlation functions of observables. This was pioneered by Rokhlin [START_REF] Rokhlin | On endomorphisms of compact commutative groups[END_REF] and has since been the focus of considerable activity. Let T be the measure preserving transformation of a dynamical system Γ on a phase space Σ. The generalized 3-strong mixing condition for any three observables or measurable sets A, B , and C take the form lim m,n→∞

µ A ∩ T m B ∩ T m+n C = µ (A) µ (B) µ (C) (4.1)
The strong mixing condition for the more general case of l-fold mixing takes the form lim

α 1 ,...,α k →∞ µ A 1 ∩ T α 1 A 2 . . . T k-1 i=1 α i A i = k i=1 µ (A i ) (4.2)
In the equivalent language of correlations functions we define the l-fold mixing correlation function of l + 1 observables, f i (x), i = 1, 2, . . . , l + 1, l = 1, 2, . . . , x ∈ Σ as follows:

C n 1 ,n 2 ,...,n l (f 1 , . . . , f l+1 ) = Σ dµ(x)f 1 (x)f 2 (T n 1 x)f 3 (T n 1 +n 2 x) • • • f l+1 (T n 1 +n 2 +•••+n l x)- l+1 i=1 Σ dµ(x) f i (x) (4.3)
where n 1 , n 2 , . . . , n l = 1, 2, . . . We say that the dynamical system (Γ, T, Σ) exhibits l-fold mixing iff

C n 1 ,n 2 ,...,n l (f 1 , . . . , f l+1 ) → 0 (4.4) as (n 1 , n 2 , . . . , n l ) → ∞, for observables {f i (x)} i=1,...,l+1 which are smooth enough and square integrable, as well as all their derivatives. We are interested in the case when T is the Arnol'd cat map and Σ = T 2 , the two-dimensional torus with radii equal to 1 and measure dµ(x) = d 2 x. In this section we shall show that the ACM exhibits l-fold mixing for every integer l = 1, 2, . . . and compute the corresponding mixing times, generalizing the calculations of the previous section.

To this end, we expand the f i (x) in Fourier series,

f i (x) = k∈Z×Z c (i) k e 2πik•x (4.5)
and extract the constant part:

f i (x) = c (i) 0 + f i (x) (4.6)
where c

(i) 0 = T 2 d 2 x f i (x) (4.7)
and

f i (x) = k =(0,0) c (i) k e 2πik•x (4.8)
We apply this decomposition only to f 1 (x) and obtain the following expression for C n 1 ,n 2 ,...,n l (f 1 , . . . , f l+1 ) :

C n 1 ,n 2 ,...,n l (f 1 , . . . , f l+1 ) = c (1) 0 Σ d 2 x f 2 (T n 1 x) • • • f l+1 (T n 1 +•••+n l x)+ d 2 x f 1 (x)f 2 (T n 1 x) • • • f l+1 (T n 1 +•••+n l x) - l+1 i=1 c (i) 0 (4.9)
The first integral in this expression can be rewritten, upon performing the change of variables, x → y = T n 1 x, that leaves the measure invariant, as

Σ d 2 x f 2 (T n 1 x) • • • f l+1 (T n 1 +•••+n l x) = Σ d 2 x f 2 (x) • • • f l+1 (T n 2 +•••+n l x) = C n 2 ,n 3 ,...,n l (f 2 , f 3 , . . . , f l+1 ) + l+1 i=2 c (i) 0 (4.10)
whence we deduce the "recurrence relation"

C n 1 ,...,n l (f 1 , f 2 , . . . , f l+1 ) = c (1) 0 C n 2 ,n 3 ,...,n l (f 2 , f 3 , . . . , f l+1 ) + d 2 x f 1 (x)f 2 (T n 1 x) (4.11)
What is noteworthy is that the product of the constant terms has been eliminated.

If we now assume that C n 2 ,...,n l (f 2 , . . . , f l+1 )-the l -1st order correlation function, tends to 0 as n 2 , . . . , n l → ∞, it remains to prove that the last term in eq. (4.11), also, tends to 0, as n 1 , n 2 , . . . , n l → ∞.

To this end we shall employ the method of the previous section, making use of the selfadjoint differential operator,

D u = 1 2πi u • ∂ (4.12)
where ∂ = (∂ x 1 , ∂ x 2 ) and u = (u 1 , u 2 ) is the eigenvector of the ACM that corresponds to the largest eigenvalue of the ACM. The inverse differential operator, [D u ] -1 acts on the complex functions of the torus that don't have a constant term. This is the reason it's useful to extract it by writing the function f 1 (x) as f 1 (x) = c

(1) 0 + f 1 (x). We notice, now, that

I n 1 ,...,n l ≡ T 2 d 2 x f 1 (x)f 2 (T n 1 x) • • • f l+1 (T n 1 +...+n l x) = - T 2 d 2 x [D u ] -1 f 1 (x) D u (f 2 • • • f l+1 ) (4.13) where D u f i (x) = k∈Z×Z c (i) k (k • u)e 2πik•x (4.14)
for i = 2, . . . , l + 1 and

D -1 u f 1 (x) = k∈Z×Z c (1) k 1 k • u e 2πik•T n x (4.15)
as well as

D u f i (T n x) = λ n k∈Z×Z c (i) k (k • u)e 2πik•x (4.16)
We can now evaluate the action of D u on the product

f 2 • • • f l+1 ; we find D u f 2 (T n 1 +n 2 x)f 3 (T n 1 +n 2 +n 3 x) • • • f l+1 (T n 1 +n 2 +•••+n l+1 x) = l+1 m=2 F l,m (x)D u f m (T n 1 +n 2 +•••+nm x) (4.17) where F l,m (x) = l+1 i=2,i =m f i (T n 1 +n 2 +•••+n i x) (4.18)
These imply that eq. ( 4. [START_REF] Kosenko | On Rokhlin's multiple mixing problem[END_REF]) can be written as [START_REF] Kosenko | On Rokhlin's multiple mixing problem[END_REF]) In order to obtain the desired result, that the integral in eq. (4.13) does vanish-and, what's much more interesting, how does it vanish-in the limit n 2 , n 3 , . . . , n l → ∞, we apply the triangle and the Cauchy-Scharz inequalities repeatedly:

D u f 2 (T n 1 +n 2 x)f 3 (T n 1 +n 2 +n 3 x) • • • f l+1 (T n 1 +n 2 +•••+n l+1 x) = l+1 m=2 F lm (x)ρ n 1 +n 2 +•••+nm f ′ m (T n 1 +n 2 +•••+nm x) (4.
|I n 1 ,n 2 ,...,n l | 2 ≤ ||D -1 u f 1 || 2 ||D u [f 2 • • • f l+1 ]|| 2 ≤ ||D -1 u f 1 || 2 l+1 m=2 ρ 2(n 1 +•••+nm) ||F lm (x)|| 2 ||f ′ m || 2 (4.20)
The norm of F lm can be bounded as

||F lm (x)|| 2 ≤ l+1 i=2,i =m ||f • T n 1 +n 2 +•••+n i-1 x|| 2 (4.21)
We notice that due to the measure preserving properties of T the various factors do not depend on n 1 , n 2 , . . . , n l :

||f • T n 1 +n 2 +•••+n i-1 x|| 2 = d 2 x |f i (T n 1 +n 2 +•••+n i-1 x)| 2 = d 2 x |f i (x)| 2 (4.22)
The same occurs for the ||f ′ m || 2 .

Collecting together all terms we obtain the effective bound using the eigenvector u -, corresponding to the eigenvalue ρ -< 1, therefore

|I n 1 ,...,n l | 2 ≤ l+1 m=1 d m ρ 2(n 1 +n 2 +•••+nm) - (4.23)
where d m contain all numerical factors.

To leading order, therefore, we find that

|I n 1 ,...,n l | 2 → d 1 ρ 2n 1 -= d 1 e -2n 1 log λ + (4.24)
thereby completing the induction hypothesis. That this term, indeed, vanishes, in the long time limit, guarantees the "sufficient" part of Rokhlin's conjecture, that 2-mixing induces l-fold mixing for any l > 2.

The mixing time for n symplectically coupled Arnol'd cat maps

It is, now, interesting to examine the mixing properties of the system of n > 1 coupled Arnol'd cat maps. It is known that hyperbolic (Anosov) linear maps on compact phase spaces exhibit strong, as well as l-fold (for all l = 1, 2, . . .), mixing. In this section, therefore, we shall focus on the effects the coupling has on the mixing properties of such maps, generalizing the calculation for the single map, that was studied in section 3 (the case of l-fold mixing is a straightforward generalization).

The idea is to use the differential operators,

D i u = 1 2πi u i • ∂ = 1 2πi u i a ∂ a (5.1)
that are the generalization for n degrees of freedom of the differential operators (4.12) that were used for computing the mixing time of one cat map. Here ∂ a = ∂/∂x a and a = 1, 2, . . . , 2n labels the point on the 2n-dimensional torus. u i is an eigenvector of the evolution operator

M = I C C I + C 2 (5.2)
where the symmetric, n × n integer matrix C is parametrized as

C = K + G(P + P T ) (5.3)
for the case of the closed chain of n maps and K and G are diagonal matrices.

where S K-S is the Kolmogorov-Sinai entropy of the system of n coupled maps. Therefore, we conclude that |C r (f, g * )| ∼ e -rS K-S , (5.11) as r → ∞, and this allows us to identify the mixing time of the system of n coupled ACM maps with 1/S K-S :

τ mixing (ACM lattice) = 1 S K-S

(5.12)

The steps of this computation can be straightfowardly generalized to the case of so-called "l-fold mixing", as was done for the single cat map in the previous section. These bounds hold even for the case of observables with square integrable first order partial derivatives. In the case that higher order derivatives are square integrable, we can use correspondingly higher powers of the differential operators D (i) uand obtain, correspondingly, better bounds.

Conclusions and outlook

In this work we have studied, in detail, the large time asymptotic behavior of the correlation functions of observables, that describe the mixing properties of a special class of automorphisms, those of coupled Arnol'd cat maps, which are both symplectic and hyperbolic, on toroidal phase spaces and have computed the mixing time, in closed form. We have shown that the mixing time is given by 1/S K-S , where S K-S is the Kolmogorov-Sinai entropy.

In the literature numerical studies of nonlinear systems with many degrees of freedom have shown similar dependence of the mixing time on the Kolmogorov-Sinai entropy; in our case we have an analytic result, which can be used to study the problem of relaxation of chaotic systems, subject to localized initial perturbations.

At the classical level it seems that there isn't any bound on how small the mixing time of the system can be; on the other hand, at the quantum level, studies of many-body quantum systems and the evolution of localized perturbations of black hole horizons seem to indicate that such a bound does, in fact, exist and is given by log S BH , where S BH is the black hole entropy. Since the correspondance principle enforces that the short time behavior of the relaxation processes for quantum systems is given by the classical mixing properties of the system, we expect that our work is relevant for understanding the crossover of mixing from the quantum to the classical system.

Our calculation is, also, relevant for the so-called Rokhlin conjecture, for which we clarify the conditions under which it may hold.

Our results allow us to focus on the conditions for physical systems, that can minimize the mixing time, in the classical limit and set the stage for addressing how the mixing time of classical many-body systems is related to the scrambling time of the corresponding quantum systems. This is of relevance for understanding, on the one hand, transport properties of novel quantum materials, as well as the properties of quantum black holes [START_REF] Hayden | Black holes as mirrors: Quantum information in random subsystems[END_REF][START_REF] Sekino | Fast Scramblers[END_REF][START_REF] Maldacena | A bound on chaos[END_REF]. To this end it is necessary to construct the corresponding unitary evolution operators for coupled Arnol'd cat maps, going beyond our previous work [START_REF] Axenides | Modular discretization of the AdS 2 /CFT 1 holography[END_REF], for the case of the single Arnol'd cat map.

+

  and contracting along the eigenvectors, v (i) -, corresponding to the eigenvalues ρ (i) -.

  2.28)where C 0 = 0 n×n , C 1 = I n×n and C m+1 = CC m + C m-1 , with m = 1, 2, 3, . . .. This matrix recursion relation generalizes to matrices the k-Fibonacci sequence for numbers. It holds for any matrix, C and, in particular for the (symmetric, integer) matrix C, defined by eq.(2.25).The solution to this matrix recursion relation is given in terms of the Fibonacci polynomials,F m (x), with argument x = C, i.e. C m = F m (C). [13].Proof. The proof is by induction. For m = 1 it is true, by definition. If we assume it holds for m > 1, then, by the relation A m+1 = A • A m , we immediately establish that it holds for m + 1.

Acknowledgements: This research was partially supported by the CNRS International Emerging Actions program, "Chaotic behavior of closed quantum systems" under contract 318687.

As for the single map, we are interested in bounding the correlation function

for f, g ∈ L 2 (T 2n ), which have continuous and square integrable partial derivatives. We proceed as follows: We observe that to get the strong mixing property as well as the rate of convergence (the mixing time), it is enough to use a straightforward generalization of the method used in for the single map, viz. we write

(5.5) whence we, immediately, find that

Now we carry out the same procedure as for the case of one map, by acting now with the product of all the D (i)

u , i = 1, . . . , n. Upon integrating by parts as for the case of one map, we find, this time, that

(5.7)

Using the Cauchy-Schwarz inequality, we obtain the bound

(5.9)

Finally, we observe that the product over the eigenvalues can be written as


Célia Ruau

Victoria Naipal

Nathalie Gagnaire

Carlos Cantero-Martinez

Bertrand Guenet

Benoît Gabrielle
email: benoit.gabrielle@agroparistech.fr

Soil erosion has mixed effects on the environmental impacts of wheat production in a large, semi-arid Mediterranean agricultural basin

Keywords: Agricultural soils, Ebro valley, Greenhouse gas emissions, Life-cycle assessment, Soil erosion, Wheat, Erosion modelling, RUSLE, soil carbon dynamics

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Anthropogenic emissions of greenhouse gases (GHGs) need to be rapidly curbed to mitigate climate change and achieve the targets laid out in the 2015 Paris agreement [START_REF] Smith | Biophysical and economic limits to negative CO2 emissions[END_REF].

Carbon sequestration in soils has been identified as an important land based mitigation option, which is at the basis of the international "4 per 1000" initiative aiming to increase carbon storage in agricultural soils in order to partially compensate anthropogenic emissions [START_REF] Minasny | Soil carbon 4 per mille[END_REF]). It can be achieved through an increase C and nutrient inputs to soils via organic matter amendment, the establishment of cover crops or agroforestery systems, among others [START_REF] Guenet | Can N2O emissions offset the benefits from soil organic carbon storage?[END_REF]. However, a reduction in soil erosion is also necessary to preserve soil C and nutrients stocks. As a consequence, lateral fluxes due to soil erosion processes should be factored in when evaluating soil C sequestration potentials. According to the recent special report of IPCC on Climate Change and land (IPCC, 2019) a quarter of all lands are degraded and erosion is the most widespread cause for degradation. Currently, agricultural soils feature the highest soil erosion rates, which are accelerated by the removal of vegetation and ground cover, and by soil tillage practices [START_REF] Ipcc | Climate Change and Land An IPCC Special Report on climate change, desertification, land degradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems WG I WG II WG III IPCC Special Report on Climate Change, D. Lal R (2009) Soils and world food security[END_REF][START_REF] Van Oost | The impact of agricultural soil erosion on the global carbon cycle[END_REF]. Soil erosion on agricultural land has substantial impacts on soil fertility and food security (Lal 2009).

The role of erosion (mainly due to water) as a sink or source of atmospheric CO 2 is currently highly debated in the literature [START_REF] Van Oost | The impact of agricultural soil erosion on the global carbon cycle[END_REF]Lal 2009;[START_REF] Borrelli | An assessment of the global impact of 21st century land use change on soil erosion[END_REF]. Erosion causes a movement of material along the cascade of hillslopes, floodplains and rivers. This redistribution of soil material may lead, on the one hand, to C sequestration in deposition areas such as colluvial and alluvial soils due to [START_REF] Van Oost | The impact of agricultural soil erosion on the global carbon cycle[END_REF] the replacement of eroded carbon by new photosynthetic substrate on eroded land, (2) burial of eroded C in deeper soil layers, hereby reducing soil carbon respiration. On the other hand, soil decomposition can be intensified by soil erosion during sediment transport and after deposition due to aggregate breakdown, and thus lead to higher CO 2 emissions. Despite this uncertainty in erosion being a C sink or source, the above-mentioned studies showed that the effect of erosion on the C emissions cannot be ignored in C cycle research.

To estimate the impacts of erosion and its drivers, mechanistic modelling is a promising way and large efforts were done in the recent years to develop soil erosion models [START_REF] Barot | Evolving away from the linear model of research: a response to Courchamp et al[END_REF].

Nevertheless, these approaches are generally developed to represent processes occurring at a fine scale and might not be adapted to larger scales. Large scale soil erosion is usually estimated based on the Universal Soil Loss Equation Model (RUSLE) [START_REF] Panagos | A New European Slope Length and Steepness Factor (LS-Factor) for Modeling Soil Erosion by Water[END_REF][START_REF] Naipal | Global soil organic carbon removal by water erosion under climate change and land use change during AD-1850-2005[END_REF]). The RUSLE model structure also facilitates the linkage of erosion processes with the C cycle processes on land.

For example, different versions of the RUSLE model have been linked to soil organic matter models such as CENTURY (Lugato et al., 2018), and to land surface models [START_REF] Naipal | CE-DYNAM (v1), a spatially explicit, process-based carbon erosion scheme for the use in Earth system models[END_REF]. However, these modelling approached do not yet include the erosion-induced feedback on plant productivity, which may lead to an overestimation of the C input to the soil in eroding regions.

Life-cycle-assessment (LCA), is a widely used, holistic approach to estimate the environmental impacts of a given practice or product [START_REF] Hellweg | Emerging approaches, challenges and opportunities in life cycle assessment[END_REF]. It provides a suitable framework to understand the multi-faceted effects of soil erosion on GHG emissions and other environmental issues. It has been used to evaluate the impacts of land use on soil erosion [START_REF] Núñez | Inclusion of soil erosion impacts in life cycle assessment on a global scale: application to energy crops in Spain[END_REF][START_REF] Van Zelm | Spatially explicit life cycle impact assessment for soil erosion from global crop production[END_REF], or the effect of tillage practices and erosion control on the GHG balance of crops [START_REF] Martin-Gorriz | Carbon emissions and economic assessment of farm operations under different tillage practices in organic rainfed almond orchards in semiarid Mediterranean conditions[END_REF], but not on the reciprocal issue of how erosion affects the environmental performance of crops. This is understandable since most LCA of agricultural crops are done at the plot or farm scale [START_REF] Poore | Reducing food's environmental impacts through producers and consumers[END_REF], whereas to estimate the GHG balance of erosion fluxes it must be captured at the catchment scale to estimate removal and redeposition, using for instance the modelling approach described above.

Here, we applied the LCA methodology to better understand the environmental impacts of wheat production as affected by soil erosion or not at the scale of the Ebro watershed in Spain. This large catchment (85 550 km² in size) was selected because it has been exposed to erosion [START_REF] Cantero-Martínez | Long-term yield and water use efficiency under various tillage systems in Mediterranean rainfed conditions[END_REF], and cropland accounts for one third of its area. The main objective was to assess the impacts of soil erosion processes induced by rainfall and runoff on the environmental performance of the production of wheat (a common crop in the area), by combining plot-scale LCA with catchmentscale erosion modelling. The mechanistic model of soil and C erosion, CE-DYNAM, a combination of the soil erosion model RUSLE with a carbon emulator of a land surface model [START_REF] Naipal | CE-DYNAM (v1), a spatially explicit, process-based carbon erosion scheme for the use in Earth system models[END_REF] provided an estimate of net erosion (erosion losses minus deposition) from agricultural land on the hillslopes of the Ebro basin. Two crop management scenarios were compared: conservation agriculture (aimed to control erosion), referred to as CA, and conventional management (referred to as CT for conventional tillage). Yield differences were factored in using two main sources of information: plotscale observations versus a model derived from a global meta-analysis comparing CA and CT.

historically [START_REF] Cantero-Martínez | Long-term yield and water use efficiency under various tillage systems in Mediterranean rainfed conditions[END_REF][START_REF] Su | High probability of yield gain through conservation agriculture in dry regions for major staple crops[END_REF], and cropland accounts for one third of its area 2 Materials and Methods

The Ebro watershed

The watershed of the Ebro river basin extends over the autonomous communities of Cantabria, north of Castilla y Leon in the provinces of Burgos and Soria, south of Pais Vasco in the provinces of Alava, La Rioja, Navarra, Aragón, Cataluña and Castellón province in the Comunidad Valenciana.

The climate is continental-Mediterranean with hot and dry summers followed by cold winters. Main rainfall events are observed in spring (April and May) and autumn (September and October), and summer droughts occur between these periods. The Ebro basin is a closed sedimentary basin where climate, tectonics and soil erosion played an important role in its formation. Currently, the Ebro basin is characterized by tow mountain ranges, the Pyrennees in the north and Iberian System mountain ranges in the south, that surround a large alluvial plain. Therefore, soil erosion in this basin has been especially affected by the steep slopes, extensive agriculture and arid /temperate climate.

Erosion modelling

Carbon erosion and deposition flows over the Ebro watershed were estimated with the novel largescale Carbon Erosion Dynamics model Version 1 (CE-DYNAM-v1) [START_REF] Naipal | CE-DYNAM (v1), a spatially explicit, process-based carbon erosion scheme for the use in Earth system models[END_REF]. CE-DYNAM consists of a sediment budget scheme coupled to an emulator of the soil carbon scheme of the global Land Surface Model (LSM) ORCHIDEE. CE-DYNAM as a whole has been originally calibrated for the Rhine catchment, but its model structure and parameters are designed in a way to be applicable for other river catchments globally. The soil erosion module of CE-DYNAM is based on the Adjusted Revised Soil Loss Equation (RUSLE) model [START_REF] Naipal | Improving the global applicability of the RUSLE model -Adjustment of the topographical and rainfall erosivity factors[END_REF], where the yearly average soil erosion removal is calculated as a product of the topographical factor (S), the rainfall erosivity factor (R), the land cover and management factor (Cm), and the soil erodibility factor (K).

The Cm factor includes the effects of vegetation cover and crop residues on soil erosion for tilled soil (conventional agriculture). The effect of conservation practices on erosion are not taken into account here. The simulated yearly soil erosion rate (E) represents the potential gross removal of soil from eroding hillslopes.

The daily gross carbon removal rate by soil erosion (Ce) can then be calculated as following:

where BD top is the bulk density of the surface layer (g cm -3 ) and CER is the enrichment ratio, i.e. the ratio of the carbon content in the eroded soil to that of the source soil material, set equal to 1. dz is the soil layer thickness in meters. CE-DYNAM simulates the dynamics of soil organic C (SOC) down to the 2 m depth.

A fraction (f) of the eroded soil gets deposited at the foot of the hillslopes (colluvial deposits). This fraction is a function of the slope gradient and the vegetation type.

Where a f and b f are constant parameters that relate f to the average topographical slope (θ) of a grid cell depending on the type of land cover. These parameters have been derived from the study of [START_REF] Naipal | CE-DYNAM (v1), a spatially explicit, process-based carbon erosion scheme for the use in Earth system models[END_REF]. θ max is the maximum topographical slope of the watershed. The rest of the eroded flux is transported to the floodplains and water network. In this study we will only consider the soil erosion and deposition, with resulting C dynamics on hillslopes. CE-DYNAM has three SOC pools, labelled active, slow and passive, respectively. These pools differ by their C residence time, which is shortest for the active pool and longest for the passive one. Changes in each SOC pool on the hillslopes are a function of the C input to the soil by litter decomposition (I), the C decomposition rate (k), which is dependent on the soil physical parameters such as temperature, and the net erosion rate (k E ), which is 1 minus the gross soil erosion (E) times the deposition rate (f). The I and k are landcover type specific and derived from the ORCHIDEE LSM on which forms the basis of the SOC scheme of CE-DYNAM.

For a more detailed description of the model see [START_REF] Naipal | CE-DYNAM (v1), a spatially explicit, process-based carbon erosion scheme for the use in Earth system models[END_REF]. between the erosion and no-erosion simulations. Here we consider the C fluxes provided by CE-DYNAM when erosion is activated as a surrogate of conventional agriculture practices whereas simulations performed when erosion is not activated are surrogate for conservation agriculture practices that limit erosion fluxes.

Input data for CE-DYNAM

The input datasets for CE-DYNAM as described below are the same as in [START_REF] Naipal | CE-DYNAM (v1), a spatially explicit, process-based carbon erosion scheme for the use in Earth system models[END_REF].

The land cover fractions are derived from the historical 0.25• maps of [START_REF] Peng | Sensitivity of land use change emission estimates to historical land use and land cover mapping[END_REF]. Daily precipitation data for the period 1850-2005 is derived from the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP), product ISIMIP2b Frieler et al., (2017). These data are based on model output of the Coupled Model Intercomparison Project. Phase 5 (CMIP5 output of IPSL-CM5A-LR [START_REF] Taylor | An overview of CMIP5 and the experiment design[END_REF], which are bias corrected using observational datasets and the method of [START_REF] Hempel | A trend-preserving bias correction : The ISI-MIP approach[END_REF] and made available at a resolution of 0.5 •. Data on soil bulk density and other soil parameters to calculate the soil erodibility factor (K), available at the resolution of 1 km, have been taken from the Global Soil Dataset for use in Earth System Models (GSDE) [START_REF] Shangguan | A global soil data set for earth system modeling[END_REF]. Finally, the slope steepness factor (S), was originally estimated at the resolution of 1 km, based on the GTOPO30 digital elevation model. All datasets that were not available at a 5 arcminute resolution were regridded to this specific spatial model resolution before calculating the erosion rates.

Life-Cycle Assessment

Goals and scope

The objectives of the LCA were two-fold: to appraise the contribution of soil erosion processes to the GHG balance of wheat, and to examine the effect of erosion control practices on the environmental performance of this crop in the context of the Ebro basin. To single out the effect of soil erosion, two sets of simulations were carried out with the catchment-scale model CE-DYNAM: a scenario with minimum erosion (reflecting erosion control practices) and a scenario with standard, default erosion (reflecting conventional crop management practices thus far). These simulations correspond to "eroded" and "non-eroded" watersheds, respectively, associated with CT and CA management practices.

To compensate for the loss of crop productivity associated with the losses of soil material and nutrients, fertilizer input rates were adjusted in the eroded case based on the estimated rates of N and P losses. These were estimated from the net C losses output by the erosion model assuming C:N and C:P ratios of 10 and 38, respectively, in eroded material [START_REF] Xu | A global analysis of soil microbial biomass carbon, nitrogen and phosphorus in terrestrial ecosystems[END_REF]. Nevertheless, we did not represent the effect of soil erosion on physical and hydraulic parameters (e.g. porosity, water holding capacity).

We used cradle-to-farm-gate system boundaries and the functional unit for this system was the production of one kilogram of wheat. The time-frame considered in the inventory extended from the harvest of the previous crop to the harvest of wheat.

Life-cycle inventory (LCI)

LCA was carried out with the SimaPro software (SimaPro 8.0.3, Pré Sustainability, LE Amersfoort, NL), and two LCI databases: AGRIBALYSE [START_REF] Colomb | Agribalyse, the French LCI database for agricultural products: high quality data for producers and environmental labelling[END_REF], which focuses on agricultural products in France, and EcoInvent (v3.6; Ecoinvent, Zurich) for background processes.

Wheat production was adapted to the Ebro valley conditions based on the management and yield data of [START_REF] Cantero-Martínez | Long-term yield and water use efficiency under various tillage systems in Mediterranean rainfed conditions[END_REF], which reports on various tillage systems at 3 field sites in the Ebro river valley. The no-till system was used for the non-eroded watershed modelling while a cropping system involving mouldboard plough was associated with the eroded watershed. The "Agramunt" site was chosen since it involved wheat crops every other year from 2000 to 2005. Crops received organic fertilizer inputs in this trial. However, since inorganic fertilizer types predominate in the Ebro catchment, organic input rates were converted into equivalent mineral fertiliser inputs based on the fertilising values of manure and slurry (Table 1).

Wheat yields are subject to large inter-annual and spatial variations across the Ebro basin (Cantero-Martinez et al., 2007;[START_REF] Plaza-Bonilla | No-tillage reduces long-term yield-scaled soil nitrous oxide emissions in rainfed mediterranean agroecosystems: A field and modelling approach[END_REF]. Since yield plays an important role in the outcome of LCAs, we used various sources of information and scenarios to determine crop yields for the conservation and standard tillage practices.

In the first yield scenario, the wheat yield in conventional agriculture (eroded watershed) was calculated from the average of the wheat yields of each province and autonomous community in the Ebro catchment area in Spain from 1996 to 2006. Yield data were taken from the annual report of the Ministry of Agriculture, Fisheries and Food of the Spanish Government. The average wheat yield in the eroded watershed analysis was 3283 kg grains ha _1 . Then, to assess the yield in conservation agriculture, we used the results of a model built from a meta-analysis comparing the yields of conventional (CT) and conservation agriculture (CA) wheat worldwide [START_REF] Su | High probability of yield gain through conservation agriculture in dry regions for major staple crops[END_REF], using machine-learning and a global data-base. The average yield change ratio between CA and CT from 1996 to 2006 was estimated at -0.010 over the Ebro watershed area. This led to an estimate of 3250 kg grains ha -1 for the yield of wheat under conservation agriculture practices. This scenario has the advantage of relying on catchment-scale data and a global model. Nevertheless, yields vary widely depending on climate, soil conditions and crop management characteristics. In some particular cases, conservation agriculture achieves higher yields than its conventional counterpart, as reported by [START_REF] Lampurlanés | Long-term analysis of soil water conservation and crop yield under different tillage systems in Mediterranean rainfed conditions[END_REF] in the Ebro basin. We therefore developed a second yield scenario based on these observations One of the main reasons for this is that CA not only reduces erosion but also improves the water-holding capacity of the soil. The average wheat yield considered in this scenario was averaged over 7 years at the studied region and was 2314 kg ha -1 for the no-till system (non-eroded watershed) and 1436 kg ha -1 for the CT system (eroded watershed). [START_REF] Lampurlanés | Long-term analysis of soil water conservation and crop yield under different tillage systems in Mediterranean rainfed conditions[END_REF] in the Ebro basin. We therefore developed a second yield scenario based on these observations One of the main reasons for this is that CA not only reduces erosion but also improves the water-holding capacity of the soil. The average wheat yield considered in this scenario was averaged over 7 years at the studied region and was 2314 kg ha -1 for the no-till system (non-eroded watershed) and 1436 kg ha -1 for the CT system (eroded watershed).

The net effect of erosion on soil C variations at the field scale was derived from the transient CE-DYNAM simulations, averaged over the 1985-2005 time slice. This erosion-induced C sink/source term for hillslopes compounded the following processes and fluxes: C replacement on eroding soils by new litter input (photosynthesis), C burial in deposition sites, increased decomposition of deposited SOC and the lateral C flux that leaves the hillslopes of the Ebro basin. The floodplains were excluded from the analysis, and it was assumed that the C flux leaving the hillslopes was entirely respired and lost.

Impact characterization

Environmental impact assessment was carried out with the IMPACT World+ characterization method (IMPact Assessment of Chemical Toxics), as integrated in the SimaPro software [START_REF] Bulle | Impact World+: a globally regionalized life cycle impact assessment method[END_REF]. This method allows analysing the results with 15 midpoint categories which are linked to 4 damage categories. The following midpoint categories were included: human toxicity (cancer and noncancer), respiratory inorganics, ionizing radiation, ozone layer depletion, respiratory organics, aquatic ecotoxicity, terrestrial acidification, land occupation, aquatic acidification, aquatic eutrophication, global warming, non-renewable energy and mineral extraction.

3 Results and discussion

Evaluation of the erosion model

We show that erosion has an important impact on the total SOC stock averaged for the period 1995-2005 in crop fields at 30cm depth (Fig. 2a) with a total SOC stock of 86.5 t C ha -1 when erosion was set to zero and 11.4 t C ha -1 when erosion was considered. It must be noted here that the big differences between the two simulations is partially due to the fact that maximized erosion is considered during the entire simulations (1860-2005) whereas it was calibrated on modern agriculture techniques. Thus, we probably overestimate the erosion fluxes for the first decades of the simulation since mechanization of agriculture responsible for higher erosion rates started in this region in the 1940's. To evaluate the model we extracted data on SOC stocks from the GSDE product [START_REF] Shangguan | A global soil data set for earth system modeling[END_REF]) when plant cover was at least 30% crop using the land cover maps from [START_REF] Peng | Sensitivity of land use change emission estimates to historical land use and land cover mapping[END_REF]. We calculated over the studied crop regions an average SOC stock at 30cm depth of 67.3 t C ha -1 . This lies in between the values predicted by the model with or without erosion, knowing that the model simulations are extreme cases considering no erosion at all or erosion rates calibrated over conventional agriculture without taking into account the erosion control measures that can be locally deployed in the fields. We also compared our simulated erosion rates over the Ebro basin with findings of other studies [START_REF] Quine | Soil erosion and redistribution on cultivated and uncultivated land near las bardenas in the central Ebro river Basin, Spain[END_REF][START_REF] Lizaga | Estimating soil redistribution patterns with 137Cs measurements in a Mediterranean mountain catchment affected by land abandonment[END_REF]). In Fig. 2b explained the differences between the two studied regions. However, also the fact that we assume a model steady state in the period 1850-1860 for CE-DYNAM may contribute to the underestimation in erosion. It is important to take into account the past changes in climate and land use as possible in order to be able to better predict the current state of soil erosion but the absence of precise knowledge on previous land use and land management makes this task very complex. Furthermore, the SOC is derived from a simulation performed by the coarse resolution ORCHIDEE LSM and the plant cover is derived from a global map, which may also contribute to the discrepancies. Finally, net soil erosion rate from CE-DYNAM can also be compared to observations derived from [START_REF] Lizaga | Estimating soil redistribution patterns with 137Cs measurements in a Mediterranean mountain catchment affected by land abandonment[END_REF] and from [START_REF] Quine | Soil erosion and redistribution on cultivated and uncultivated land near las bardenas in the central Ebro river Basin, Spain[END_REF]. In this case the model predicts a net soil erosion rate ranging from 1.7 to 4.7 t ha -1 yr -1 (equivalent to a range from 43 to 95.5 kg C ha -1 yr -1 ) which is lower than the rate derived from [START_REF] Lizaga | Estimating soil redistribution patterns with 137Cs measurements in a Mediterranean mountain catchment affected by land abandonment[END_REF] who estimated a net soil erosion rate of 8.3 t ha -1 yr -1 . This could be expected as the gross erosion rate is already underestimated by the CE-DYNAM model. [START_REF] Quine | Soil erosion and redistribution on cultivated and uncultivated land near las bardenas in the central Ebro river Basin, Spain[END_REF], also using the 137 Cs tracer method, estimated a rate ranging from 16 to 25 t ha -1 yr -1 for a small sub basin of 0.52 km 2 . The representativeness of such sub basins compared to the entire Ebro basin is therefore not guaranteed. Finally, since conservation agriculture practices are not implemented in the CE-DYNAM model that can run with no erosion or with maximized erosion for every grid cells, the performances of the model are reasonable. Net erosion rates might be under estimated by CE-DYNAM but the absence of data over the entire Ebro basin makes the evaluation of the model not straightforward.

Contribution of soil erosion to the GHG balance of wheat

Based on the CE-DYNAM simulations of soil organic C variations over the recent time slice, cropland soils emerged as net sources of CO 2 , with an estimated magnitude of 215 and 249 kg CO 2 ha -1 yr -1 for the non-eroded and eroded watersheds, respectively. It is notable that despite the large effect of erosion on soil C stocks (Fig. 2), the losses experienced since the 1860 starting point in the eroded watershed were still more than compensated by the rates of C burial in deposition sites.

Indeed, erosion drastically affect the spatial distribution of SOC with higher SOC stocks in floodplain and lower SOC stocks in hillslopes. Ultimately, this sink effect was no longer active in the non-eroded configuration.

Figure 3 puts these soil C sources in perspective by reporting the life-cycle GHG emissions of wheat production for yield scenario #1, per kg of wheat grains harvested. Overall, they only accounted for 5% of life-cycle GHG emissions. Their magnitude fits with the 0 to 1700 kg CO 2 ha -1 yr -1 range reported by [START_REF] Goglio | A comparison of methods to quantify greenhouse gas 1 6 emissions of cropping systems in LCA[END_REF] in their survey of the life-cycle GHG balances of wheat crops worldwide, although some studies based on direct field observations of SOC variations pointed to larger sink or source strengths (from -3.5 to 5 Mg CO 2 ha -1 yr -1 ). Conservation agriculture in general is assumed to favour soil C sequestration, and Alvaro-Fuentes et al. ( 2011) estimated it could store 1 1

between 200 and 500 kg C ha -1 yr -1 (730 to 1 830 kg CO 2 ha -1 yr -1 ) whereas conventional systems were simply neutral. The fact that erosion simulations started prior to 1985 (the first year of the 20-year sequence used here for SOC variation calculations) and that CA and CT systems had largely different SOC contents at that time (reflecting the "legacy" of past erosion processes) certainly mitigated the sink capacity of CA. Factoring in the past (1860 to 1985) differences in SOC simulated with CE-DYNAM between the two systems would amount to attributing an extra 480 kg C ha -1 yr -1 sink term to the CA system. This would offset 40% of the life-cycle GHG emissions of the CA systems, which are equivalent to 1210 kg C ha -1 yr -1 . Along with potential soil C sink terms, Figure 4 highlights the importance of N 2 O emissions (which made up 85% of direct field emissions of GHG), and the use of chemical fertilizers on the GHG balance of wheat. Nitrogen fertilizers bore the largest impact, followed by phosphorus and potassium fertilizers. The 35% share of these inputs in the overall carbon footprint was similar to the 30% ratio reported by [START_REF] Poore | Reducing food's environmental impacts through producers and consumers[END_REF]. From a fertilizer use perspective, another consequence of erosion is the need to compensate for nutrient losses. This adjustment has a two-fold impact: a rise in off-farm emissions due to the manufacturing and transportation of extra inputs, and an increase in the field emissions of N 2 O derived from fertilizer N.

Compared to the non-eroded scenario, the eroded case resulted in a 5% relative increase for both terms, and thus a 10% relative increase overall (Figure 3).

Wheat production with CT practices required more frequent use of agricultural machinery (ploughing, fertilization) inducing a slight increase in CO 2 emissions. Indeed, diesel consumption was 3-times higher with CT than with CA, but contributed only a minor fraction (under 10%) of the emissions. Also, erosion per se lead to higher impacts for some of the impact categories due to the pollution of surface water through runoff of nutrients, organic matter and chemicals. This was particularly apparent for freshwater eutrophication and ecotoxicity, which were 25% lower under CA (Figure 3 and Table 2). Conversely, the loss of soil with CT mitigated to a minor extent (5%) the impact of CA wheat on human toxicity (non cancer) because a reduced exposure of human consumers to contaminated soils.

Here, we considered that only mineral forms of fertilizers inputs whereas organic ones may impact the productivity and the SOC dynamics of cropland. For instance, [START_REF] Larney | A simulated erosion approach to assess rates of cattle manure and phosphorus fertilizer for restoring productivity to eroded soils[END_REF] showed that organic fertilizers were more effective at supplying P to crops in eroded areas than the input of phosphates, which were immobilized by calcium carbonates. Comparing inorganic and organic fertilizers to maintain crop productivity on eroded lands could bring interesting insight into the environmental performance of cereals in eroded catchments.

Erosion and crop yields

Figure 3 compares the environmental impacts of wheat production under CA and CT using a first scenario in which yields differed by only 1% between both managements. Under this assumption of 1 2 similar yields, soil erosion per se caused higher emissions to the air and water compartments in general, and resulted in 10 to 25% larger impacts for all categories except human toxicity. With the second scenario based on cropping systems trials comparing CA and CT, differences between the two managements and their associated erosion configurations were much larger, reaching 60 to 80% (Fig. 3b), again with the exception of human toxicity (non cancer). This magnitude could be expected from the 61% relative yield difference between CA and CT wheat (Table 1), and the differences in environmental impacts reported in the previous sections, whereby CA outperformed CT because of soil erosion processes.

Overall, experimental trials on the effect of soil erosion rates on crop yields are inconsistent, but a recent meta-analysis pointed out that they were negative when erosion depth dropped under 10 cm [START_REF] Zhang | Effect of soil erosion depth on crop yield based on topsoil removal method: a meta-analysis[END_REF]). In the Ebro valley, the soil losses simulated by CE-DYNAM in the eroded watershed (averaging 1.5 t ha -1 yr -1 ) indicate a rapid depletion of the topsoil and point at an erosion depth larger than this threshold (see also picture on Figure 1). This adverse effect on crop productivity is backed by the field data used for yield scenario 2, with the additional benefits of soil conservation practices for the CA/non eroded watershed configuration. These benefits are known to be larger under semi-arid conditions [START_REF] Su | High probability of yield gain through conservation agriculture in dry regions for major staple crops[END_REF], even though the model we used only pointed to a marginal effect. The CA system also featured less intensive management practices than its conventional counterpart, related to the absence of tillage and reduced fertilizer inputs. No-till is known to trade-off with other practices such as weed management, and collecting management data representative of the area is paramount to produce reliable results.

Finally, an important point to mention is that our estimation of erosion from CE-DYNAM does not consider the feedbacks between erosion and primary productivity. CE-DYNAM is a C-only model and the loss of soil due to erosion only remove C whereas in reality soil erosion also remove N and P and affect soil physical structure, water infiltration etc. Therefore, the higher SOC loss due to erosion over the 1985-2005 period given by CE-DYNAM can be underestimated because primary production and therefore C inputs into the soil is not affected by erosion. Our estimation of erosion effect on SOC should be therefore considered as an optimistic estimation likely underestimated.

The full GHG balance of erosion

Here we demonstrated that soil erosion increased the GHG emissions of wheat through its impact on SOC dynamics, nutrient losses, chemical fertilizer use and crop management practices and, most importantly, crop yields. The difference ranged between 0.1 and 2.5 kg eq. CO 2 kg -1 of wheat grains produced depending on the yield scenario. This is quite significant considering the world average of 1.1 kg eq. CO 2 kg -1 reported by [START_REF] Poore | Reducing food's environmental impacts through producers and consumers[END_REF] for this crop, and the average value of 0.89 1.1 kg eq. CO 2 kg -1 reported by the AgriBalyse database in France [START_REF] Colomb | Agribalyse, the French LCI database for agricultural products: high quality data for producers and environmental labelling[END_REF]. Compared to 1 3

the control watershed, erosion on its own increased on-field emissions of CO 2 (related to SOC losses), indirect emissions from chemical fertilizer use while diesel consumption was higher with conventional crop management than with soil conservation practices. This resulted in a compounded increase of 7% in the GHG emissions of wheat production, yields effects notwithstanding. On the other hand, erosion leads to the burying of soil organic carbon and represents a net sink [START_REF] Van Oost | The impact of agricultural soil erosion on the global carbon cycle[END_REF]. Using the simulations done by CE-DYNAM we calculated this sink to be 6.7 g C m -2 yr -1 , which corresponds to 5.5% of the wheat's life-cycle emissions. Factoring in this sink to obtain the "full" GHG balance of the eroded scenario, we find that the effect of C burying soil is commensurate with the GHG emissions due to other erosion-related processes (which a share of 7%).

Whether erosion is a net sink or a source of C is still debated (see [START_REF] Doetterl | Erosion, deposition and soil carbon: A review ofprocesslevel controls, experimental tools and models to address C cycling in dynamic landscapes[END_REF] for review) but its effects on the life-cycle impacts of crops had hitherto not been taken into account. We show that in the Ebro watershed, even though a fraction of soil C ends up buried in river sediments through erosion, this fraction is on a part with the GHG emissions directly caused by erosion, and small compared to the overall life-cycle emissions of wheat, making erosion carbon-neutral at best.

Conclusion

Soil erosion had mixed effects on agricultural GHG emissions over the Ebro watershed. From a life-cycle perspective, taking into account all the processes involved in wheat production in this area, erosion was both a source of emissions and a sink through the burial of C in river sediments. These two terms being of similar magnitude, erosion per se appeared as a carbon-neutral process. The major differences between the eroded and non-eroded watersheds occurred with the scenario assuming a large yield gap between the conventional and conservation-orientated crop management practices. In this case the zero-erosion control clearly outperformed its eroded counterpart, with 60 to 80% lower life-cycle impacts per kg of wheat produced. Notwithstanding yield differences, soils in the eroded watershed had lost a large amount of carbon by the beginning of the time period investigated here (1985 to 2005), leading to an overall loss of c. 0.5 t C ha -1 yr -1 relative to the un-eroded control, potentially offsetting 40% of the wheat's C footprint. Factoring in this legacy of past land use changes and management practices is not an easy task in retrospect. However, our modelling and assessment

shows that it is important to control erosion and to adopt conservation practices to preserve soil structure, crop productivity, and mitigate the environmental impacts of crops in the context of climate change.

  CE-DYNAM was run for the period 1850-2005, with and without soil erosion dynamics, with a daily time-step at a spatial resolution of 5 arcminutes (~ 8km). First the model was run to a steady state based on the environmental conditions of the period 1850-1860, and afterwards transient simulations were performed with changing climate and land use conditions. The overall change in SOC stocks over the period 1850-2005 is then a result of climatic changes, increasing atmospheric CO 2 concentrations, land use change and soil erosion. The resulting net soil C emissions due to erosion can be calculated as the difference in the cumulative SOC change over the period 1850-2005

Figure 1 :Figure 2 : 4 Figure 3 : 5 Figure 4 :

 124354 Figure 1: Evidence of soil erosion processes on agricultural land in the Ebro valley (credit: Julien Tournebize, INRAE).

  , we compared the gross erosion rate predicted by the CE-DYNAM model and the rate estimated by[START_REF] Lizaga | Estimating soil redistribution patterns with 137Cs measurements in a Mediterranean mountain catchment affected by land abandonment[END_REF] using137 Cs tracer. Our model underestimates the gross erosion rate compared to the rate derived from Lizaga

	et al., (2018). The model predicts a gross soil erosion rate ranging from 6.6 to 18.2 t ha -1 yr -1
	(equivalent to a range from 167.8 to 371.7 kg C ha -1 yr -1 ), whereas Lizaga et al., (2018) estimated a
	gross soil erosion rate over their studied region to be 25.4 t ha -1 yr -1 . It is important to note that Lizaga

et al., (2018)

did not estimate the soil gross erosion rate over the full Ebro basin region but over a sub basin of 23 km 2 located in Navarra province that is a Mediterranean-continental climate suffering of Atlantic influence receiving more rainfall that in the central Ebro valley area. This sub basin is also generally steeper than the central Ebro valley area. The model-observation errors may therefore be 1 0

1 9 Table 1 .

 91 Crop management data for wheat production on the non-eroded and eroded watersheds (corresponding to conventional tillage -CTand conservation agriculture -CApractices, respectively), for yield scenario #1. The extra N and P inputs for CT compared to CA compensate for the net loss of organic matter through erosion.

			Non-eroded watershed (CA)	Eroded watershed (CT)
			Barley from 1990 to 2000 then wheat one year in two until 2005
		Crop rotation
			3250 kg grain DM ha -1	3283 kg grain DM ha -1
		Yield in scenario
		#1
		Yield in scenario	2314 kg grain DM ha -1	1436 kg grain DM ha -1
		#2
		Tillage	No-till	1 mouldboard ploughing (30-35 cm)
				1/2 cultivator passes (15cm)
				1 cultivator pass before sowing
			No-till disc drill : 500 seeds m -2
		Sowing
			Glyphosate (0.685 kg a.i. ha -1 )	Prosulfocarb (4 kg a.i. ha -1 )
		Herbicide inputs	MCPA (0.644 kg a.i. ha -1 )
			Prosulfocarb (4 kg a.i. ha -1 )
			134 kg N ha-1	(134 + 6) kg N ha-1
		Fertilizer inputs	81.5 kg P2O5 ha-1	(81.5 + 4) kg P2O5 ha-1
		(inorganic)
			142 kg K2O ha -1	142 kg K2O ha -1
		Disease and pest	No treatment
		management
	,	Harvest	Medium-sized standard combine harvester. Straw is chopped and spread on the field.
				2
				0

Table 2 :

 2 Life-cycle impacts of the production of one kg of wheat in the non-eroded (CA) and eroded (CT) Ebro watersheds, for yield scenario #1, expressed at mid-point level.Key to units (see[START_REF] Bulle | Impact World+: a globally regionalized life cycle impact assessment method[END_REF] for more details): CTUh, CTUe: comparative toxic units for humans and ecosystems, respectively; NMVOC: non-methane volatile organic compounds.

	Impact category		Unit	Non-eroded	Eroded	watershed
				watershed (CA)	(CT)
	Human toxicity (cancer)	CTUh	2.87 10 -9	2.76 10 -9
	Human toxicity (non-	CTUh
	cancer)			3.23 10 -8	3.40 10 -8
	Particulate	matter	kg PM2.5 eq.
	formation			4.12 10 -4	4.19 10 -4
	Ionizing radiation		Bq C-14 eq.	1.8		1.9
	Ozone layer depletion	kg CFC-11 eq.	5.37 10 -8	6.08 10 -8
	Photochemical oxidants	kg NMVOC eq.
	formation			3.10 10 -3	3.63 10 -3
	Freshwater ecotoxicity	CTUe	59.2		95.8
	Terrestrial acidification	kg SO 2 eq.	1.18 10 -5	1.25 10 -5
	Freshwater acidification	kg SO 2 eq.	1.09 10 -8	1.25 10 -8
	Freshwater		kg PO 4 P-lim eq.
	eutrophication			1.13 10 -4	3.20 10 -4
	Marine eutrophication	kg N eq.	1.63 10 -4	1.74 10 -4
	Global warming		kg CO 2 eq.	1.42		1.51
	Non-renewable energy	MJ primary	5.56		6.29
	Land occupation		m2org.arable	0.206		0.204
	Mineral extraction		kg deprived resource	1.00 10 -2	1.13 10 -2
	Water scarcity		m3 world eq.	24.4		26.1

Acknowledgements

The authors would like to thank Dominique Desbois (INRAE) for supplying data on wheat yield statistics in the Ebro basin, and Yang Su (AgroParisTech/INRAE) for running his yield model for conservation agriculture on the same area. 1 5

Declarations

1 4

• Funding: this work received support from the project ERANETMED2-72-209 ASSESS funded by the European Commission. French government under the ANR "Investissements d'avenir" program with the reference CLAND ANR-16-CONV-0003 also financially supported this study.

• Conflicts of interest/Competing interests: the authors declare that they have no conflict of interest.

• Ethics approval: not applicable.

• Consent to participate: not applicable.

• Consent for publication: not applicable.

• Availability of data and material: erosion simulation data and LCA input and output data are available upon request from the authors.