Jean-Baptiste Apoung
email: jean-baptiste.apoung@math.u-psud.fr

Pascal Have
email: pascal.have@ifp.fr

Jean Houot
email: houot@iecn.u-nancy.fr

Michel Kern
email: michel.kern@inria.fr

Adrien Semin
email: adrien.semin@inria.fr

REACTIVE TRANSPORT IN POROUS MEDIA *

We present a numerical method for coupling transport with chemistry in porous media.

Our method is based on a fixed-point algorithm that enables us to coupled different transport and chemistry modules. We present the methods for solving the sub-problems, detail the formulation for the coupled problem and show numerical examples to validate the method.

Résumé. Nous présentons une méthode numérique pour le couplage du transport et de la chimie en milieu poreux. Notre méthode utilise un algorithme de point fixe, qui nous permet de coupler des modules de transport et de chimie différents. Nous présentons les méthodes numériques utilisées pour chacun des sous-problèmes, ainsi qu'une formulation pour le problème couplé, et nous validons la méthode sur quelques exemples numériques.

Introduction

Most of our drinking water comes from underground resources, but these resources are subject to numerous threats of various origins: industrial and agricultural waste, acid rains, nuclear waste, ... Understanding these risks to mitigate them, and implement remediation policies if possible, is becoming an increasingly important stake [START_REF] Steefel | Reactive transport modeling: An essential tool and a new research approach for the Earth sciences[END_REF]. In this respect, chemical phenomena play a fundamental role: chemical species in the subsurface may react both between themselves,to form new compound species, and also with the host rock. These heterogeneous reactions could potentially have an important effect on the flow and transport phenomena, by opening new channels, or closing existing ones. The mathematical model underlying these phenomena couples a transport model (advection-diffusion partial differential equations) together with a chemical model (ordinary differential equations or algebraic equations, depending on whether chemistry is modeled with kinetic or equilibrium reactions). This work is concerned with a numerical method for the simulation of coupled transport-chemistry models. We consider a simplified setting, with only equilibrium reactions, and we only take into account heterogeneous reactions in the form of sorption. The transport-chemistry problem thus takes the form of a large system of PDEs coupled with algebraic equations. As the numerical methods for both parts of the system are quite different, we want the overall coupling algorithm to respect this separation (see [START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF] for a discussion of various approaches to reactive transport simulation). In this work we use a simple fixed-point algorithm [START_REF] Mosé | Efficiency of operator splitting procedures for solving reactive transport equation[END_REF][START_REF] Montarnal | Coupling methodology within the software platform alliances[END_REF][START_REF] Saaltink | On the behavior of approaches to simulate reactive transport[END_REF], that enables us to put together the pieces for a general coupling strategy, and to demonstrate some of the phenomena involved.

An outline of the paper is as follows: in section 1 we present the models used for the basic phenomena of flow transport and chemistry. The solution method, both for the individual sub-problems and for the overall coupled problem are presented in section 2, and we show numerical examples in section 3.

Presentation of the flow, chemical phenomena and transport

We consider a porous medium occupying a domain Ω ⊂ R d (d = 2 or 3) and carrying a set of chemical species that can either be in solution (aqueous species) or attached to the host rock matrix (sorbed species). The 3 main phenomena we describe are the flow of water in the subsurface (we assume the flow is stationary, and that there is no feedback between chemistry and flow), the transport of the various species, due to this flow, and chemical phenomena both between species in solution, and between the species and the host matrix.

Flow phenomena

We consider Darcy's law and the incompressibility law to describe the velocity u and the pressure p of the flux in the domain Ω :

{ ωu = -K∇p ∇ • (ωu) = 0 in Ω, (1)
where ω is the porosity of the medium (expressing the amount of void in the porous medium), and K is the permeability tensor. The boundary conditions may be of different type

{ p = p 0 on Γ D , u.n = g on Γ N , , (2)
with p 0 a given boundary pressure, and g a prescribed boundary flux. We assume that the Darcy velocity is not affected by the chemical phenomena (i.e. flow is constant over time, and the porosity does not depend on the concentrations of the various species). We make the usual assumptions on the permeability K, which may be a tensor, and may depend on x : we assume that K ∈ L ∞ (Ω) d×d and that there exists constants k * and k * such that,

∀ξ ∈ R d , 0 < k * |ξ| 2 ≤ ξ T Kξ ≤ k * |ξ| 2 < ∞.
We also assume that the functions p 0 and g are "smooth enough".

Chemical reactions

We assume that the chemical reactions are local (that is, we consider that the reaction phenomena at a point x ∈ Ω only involve the concentrations of the species at this point). According to Rubin [START_REF] Rubin | Transport of reacting solutes in porous media: relation between mathematical nature of problem formulation and chemical nature of reaction[END_REF], chemical reactions can be classified according as to whether they are kinetic or at equilibrium, and whether they are homogeneous or heterogeneous. In this work, we only consider instantaneous reactions (there is no kinetic term), and only some forms of heterogeneous reactions. The restriction to equilibrium only reactions deserves some justification. This is a common modelling hypothesis in reactive transport in porous media (see for instance the discussions in Valocchi [START_REF] Valocchi | Validity of the local equilibrium assumption for modeling sorbing solute transport through homogeneous soils[END_REF] or Rubin [START_REF] Rubin | Transport of reacting solutes in porous media: relation between mathematical nature of problem formulation and chemical nature of reaction[END_REF]). Its validity stems from the fact that chemical phenomena usually occur on a much faster time scale than the transport phenomena, thus leading to the so-called "local equilibrium" assumption (this is of course in sharp contrast with the situation in atmospheric chemistry where taking into account kinetic reactions is essential, see for instance [3]). For the kind of applications discussed in this paper, kinetic models would be needed if we had included minerals, or biodegradation reactions. Finally, we briefly indicate at the end of section 1.4 the modifications that would be necessary if we took kinetic reactions into account.

As to the nature of the reactions, it is customary to distinguish between two types of chemical reactions:

Homogeneous reactions: These occur only within the aqueous phase, and include for instance, complexation, acid-base reactions, and oxydo-reduction reactions; Heterogeneous reactions: These obviously play a fundamental role in geochemistry. They include precipitation dissolution reactions, whereby an aqueous species can form a solid complex. These are threshold reactions, with specific difficulties, and we do not consider them further in this paper. The other type of heterogeneous reactions are the sorption reactions, which are the main vehicle of exchange at the rock-water interface. We only consider a ion-exchange mechanism, where ions are exchanged between the aqueous species and the rock matrix, and a complex is formed on the matrix. Only a specific number of "sites" on the matrix can accept ions. The medium is characterized by its Cationic Exchange Capacity (CEC). More details on the actual mechanisms involved can be found in the reference book [START_REF] Dzombak | Surface complexation modeling-Hydrous ferric oxide[END_REF]. From the simplified point of view adopted in this paper, the end result is that we can model sorption reactions in the same way as homogeneous ones. Note that the equilibrium assumption is consistent with the restriction to sorption reactions only.

We now consider N A aqueous species and N S sorbed species, undergoing N R linearly independent chemical reactions that can be written as

NA+NS ∑ j=1 A ij E j ⇆ 0, i = 1, . . . , N R .
We can express N R species, that we call secondary species, according to the Morel formalism (see [START_REF] Morel | Principes and Applications of Aquatic Chemistry[END_REF]), as functions of the N P = N A + N S -N R other species, that we call primary species.

In the following, we denote by X 1 , . . . , X Nx the concentrations of the primary aqueous species and by Y 1 , . . . , Y Ny the concentrations of the primary sorbed species (note that N x + N y = N P ). In the same way, we denote by C 1 , . . . , C Mc the concentrations of the secondary aqueous species and by S 1 , . . . , S Ms the concentrations of the secondary sorbed species (note also that M c + M s = N R ). Of course these concentrations are functions defined on R t × Ω.

To describe the chemical problem, we have two types of laws to consider :

Mass Action Law : Each reaction gives rise to a mass action law, describing how we can get the concentrations of the secondary species, given the concentrations of the primary species.

Mass action law for aqueous species : since we consider no precipitation phenomena, the i th secondary aqueous specie is obtained from a linear combination of the N x primary aqueous species. Denote by A 1,i,j the coefficient of the j th primary aqueous specie in this linear combination. The mass action law can then be written, in logarithmic form, as

log(C i ) = Nx ∑ j=1 A 1,i,j log(X j ) + log(K Ci ) (3)
where K Ci is a thermodynamic constant which depends only on the temperature (supposed constant in our model).

Mass action law for sorbed species : in the same way, the i th secondary sorbed specie is obtained from a linear combination of the N x primary aqueous species (with linear coefficients A 3,i,j ) and the N y primary sorbed species (with linear coefficients A 4,i,j ). The mass action law can be written as

log(S i ) = Nx ∑ j=1 A 3,i,j log(X j ) + Ny ∑ j=1 A 4,i,j log(Y j ) + log(K Si ) (4)
where K Si is a thermodynamic constant which depends only on the temperature.

Conservation Law : the Conservation Law (also known as the Lavoisier law) can be expressed by the two following relations

T i = Mc ∑ j=1 A 1,j,i C j + Ms ∑ j=1 A 3,j,i S j + X i , ∀i ∈ {1, . . . , N x } ( 5
)
W i = Ms ∑ j=1 A 4,j,i S j + Y i , ∀i ∈ {1, . . . , N y } (6)
where T i is the total concentration of the aqueous specie X i and W i is the total concentration (related to the CEC of the medium, and assumed constant over time in our model) of the sorbed specie Y i .

Summary : if we denote A k the matrix with coefficients A k,i,j , the Mass Action Laws (3,[START_REF] Brassard | A feasible set for chemical speciation problems[END_REF] and the Conservation Law [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Carrayrou | presentation of the MoMaS reactive transport benchmark[END_REF] can be put in the following matrix form

[ log(C) log(S) ] = [ A 1 0 A 3 A 4 ] [ log(X) log(Y ) ] + [ log(K C ) log(K S )
] and

[ T W ] = [ A T 1 A T 3 0 A T 4 ] [ C S ] + [ X Y ] (7)
where log(Z) is the vector whose i th component is equal to log(Z i ).

We can summarize these matrix equations by giving the Morel tableau:

X Y K C A 1 0 K C S A 3 A 4 K S T W (8)
This table can be easily found in the literature for concrete examples (see for example the Morel tableau given by ( 24)).

In the following, for each aqueous species i, we denote by Q i = ∑ A 1,j,i C j + X i its total aqueous concentration, and by F i = ∑ A 3,j,i S j its total fixed concentration. The associated matrix forms are

Q = A T 1 C + X and F = A T 3 S. Of course one has T = Q + F . Remark 1.1.
If the secondary species are only sorbed species (i.e. M c = 0), one can easily see that the total aqueous concentration Q i of each specie i is equal to its concentration X i .

Transport model

Since we study the transport of small particles, the mathematical model is the advection-dispersion equation (for any specie X i , Y i , C i and / or S i ) :

ω dZ dt + ∇ • (ωuZ -D(u)∇Z) = f (Z) in Ω (for aqueous species), (9)
or (as only aqueous species are subject to transport)

ω dZ dt = f (Z)
in Ω (for fixed species) [START_REF] De Marsily | Hydrogéologie quantitative[END_REF] where D(u) is a dispersion tensor given by

D(u) = ωd e I + α T ω|u| + (α L -α T )ω(u ⊗ u)|u| -1
(the Scheidegger dispersion model, see for instance [START_REF] De Marsily | Hydrogéologie quantitative[END_REF] or [START_REF] Bear | Modeling Groundwater Flow and Pollution[END_REF]), and f is a source-sink term depending on all the present species.

The scalar d e is called the molecular diffusion coefficient, α T is the transverse dispersion coefficient and α L is the longitudinal dispersion coefficient (dispersion is a macroscopic way of taking into account microscopic heterogeneities of the velocity). In practice, we often consider either molecular diffusion only or dispersion only.

The boundary conditions can be either of Dirichlet (given concentration), total flux or diffusive flux given. In the sequel, we denote the transport operator by

L(Z) = ∇ • (ωuZ -D(u)∇Z).

The coupled system

We make the classical assumption that the diffusion (and dispersivity) coefficients are the same for all species, so that the transport operator is the same. Its validity is at best questionable, but without this assumption, it would be difficult to carry out the reduction to total concentrations. We write one equation like ( 9) above for each aqueous primary and secondary species, and one equation like [START_REF] De Marsily | Hydrogéologie quantitative[END_REF] for each fixed species. Each equation involves an unknown reaction term. In order to eliminate these terms, we introduce the total aqueous concentration Q and total fixed concentration F (details on the elimination procedure are given for example in [START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF]). We obtain a transport equation of the same form for the total aqueous concentration, with a source term given by the (change in) total fixed concentration.

ω dQ dt + ω dF dt + L(Q) = 0 with L(Q) = ∇ • (ωuQ -D(u)∇Q) ( 11
)
This model is consistent with the model often used for the transport of a single species undergoing sorption (see for example [START_REF] Logan | Transport modelling in hydrogeochemical systems[END_REF]).

Finally, the transport-reaction equations can be written in the form : find C, X, S, Y , Q, F and T defined on R t × Ω such that :

[ log(C) log(S) ] = A [ log(X) log(Y ) ] + [ log(K C ) log(K S ) ] [ T W ] = A T [ C S ] + [ X Y ] T = Q + F, Q = X + A T 1 C, F = A T 3 S              (12)
ω dQ dt + ω dF dt + L(Q) = 0
Initial and boundary conditions

   (13)
It has recently been shown by S. Kraütle [START_REF] Kraütle | General multi-species reactive transport problems in porous media: Efficient numerical approaches and existence of global solutions[END_REF] that this system is well posed, under quite general hypotheses on the data. Let us note that if we had included kinetic reactions in the model, equation ( 13) would have a zeroth order term, depending on the concentrations of all species, whereas equations ( 12) would remain intact (including only equilibrium species). Thus most of our work would remain valid, the main difference being that the various transport operators would no longer be uncoupled.

Solution algorithms

In the following, we will denote by ∆t the time step (taken constant for simplicity. Variable time step could, and should,easily be included in the framework), and we wish to compute an approximation of the functions F , Q, T at time t n = n∆t, that we denote by F n , Q n and T n .

Solution of the individual problems

Flow solver

The flow computation is not the main focus of this work, but the Darcy velocity is needed for the transport. In a complex geometry, it will not be possible to prescribe a consistent velocity, so we must first compute the flow, where the main unknown will be the velocity, and not the pressure.

It has long been recognized that standard Lagrange finite elements are not adequate for this task (see [START_REF] Mosé | Application of the mixed hybrid finite element approximation in a groundwater flow model : Luxury or necessity[END_REF], [START_REF] Durlofsky | Accuracy of mixed and control volume finite element approximations to Darcy velocity and related quantities[END_REF], as the computed flux may not be continuous across a face. The standard method is the mixed finite element, based for instance on the Raviart-Thomas, or Brezzi-Douglas-Marini, spaces (see [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF] or [START_REF] Roberts | Mixed and hybrid methods[END_REF] for more details about those spaces).

We briefly review the solution of the flow problem (1) by the mixed finite element method. It is based on the mixed formulation of the system: Find (p, u) ∈ L 2 (Ω) × V g such that :

       ∫ Ω K -1 u.v dx - ∫ Ω p∇.v dx = - ∫ ΓD p 0 v.n dγ, ∀v ∈ V 0 - ∫ Ω q∇.u dx = 0, ∀q ∈ L 2 (Ω), (14)
where we have defined

V g = H g N (div, Ω) = {u ∈ L 2 (Ω), ∇.u ∈ L 2 (Ω), u.n = g on Γ N }.
Existence and uniqueness for the mixed problem [START_REF] Galassi | GNU Scientific Library Reference Manual -Revised Second Edition[END_REF] follows from the LBB theory (see [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Roberts | Mixed and hybrid methods[END_REF]), thanks to the assumptions on the data made in section 1.1.

For the discretization, we introduce a regular triangulation T h of Ω, where the elements may be triangles or rectangles in 2D, and tetrahedra or hexahedra in 3D. In the sequel, we concentrate on the rectangular case, as that is what we use in the examples. In order to obtain a positive definite linear system, we will use the hybrid form of the mixed finite element method. We first define the local Raviart-Thomas space, noting that other choices exist [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Roberts | Mixed and hybrid methods[END_REF]. Denoting by K a rectangle in the triangulation, we define

RT 0 (K) = {v h = ( α + ax β + by ) , (a, b, α, β) ∈ R 4 }. ( 15
)
0 1 00 00 11 11 00 11 00 11

w 3 w 4 E 1 w 2 w 1 E 2 E 3 E 4 1 q q 2 q 3 q 4
V g h = {v h ∈ V g , ∀K ∈ T h , v h|K ∈ RT 0 (K)}, W h = {q h ∈ L 2 (Ω), q h|K P 0 (K)},
as well as a multiplier space

Λ h = ∪ K∈T h {v h .ν K , v h ∈ RT 0 (K)}, and Λ p0 h = {λ h ∈ Λ h , ∫ ΓD (λ h -p 0 )µ h = 0, ∀µ h ∈ Λ h }
where ν K is the outgoing normal to K, and the trace is taken in the weak sense.

The mixed hybrid finite element method is obtained by relaxing the inter-element continuity for the normal fluxes, and introducing a Lagrange multiplier. Thus we look for a triple

(p h , u h , λ h ) ∈ W h × V g h × Λ p0 h such that                      ∑ K∈T h (∫ K K -1 u h .v h dx - ∫ K p h ∇.v h dx + ∫ ∂K λ h v h .n dγ, ) = 0 ∀v h ∈ V h ∑ K∈T h - ∫ K q h ∇.u h dx = 0, ∀q h ∈ W h , ∑ K∈T h ∫ ∂K µ h u h .ν dγ = 0 ∀µ h ∈ Λ 0 h . ( 16
)
The main advantage of the hybrid method is that both the pressure and the fluxes can be eliminated locally (this was the motivation for relaxing the continuity of the fluxes). Then one obtains a linear system involving only the Lagrange multipliers, which can be interpreted as pressure traces on the edges. It can be shown that this linear system is symmetric and positive definite, so that efficient solution methods can be used. It is known that the method just described is convergent, ans it has been demonstrated that it can be used for porous media, and that is gives an accurate velocity field, see the above references.

Chemical solver

Due to the widely different orders of magnitude of the concentrations that are commonly encountered (this will be seen in section 3.1 below, for the MoMaS benchmark [START_REF] Momas | Exercices de qualification numérique de codes[END_REF], or the gallic acid example [START_REF] Brassard | A feasible set for chemical speciation problems[END_REF]), we solve this system by taking as unknown x = log(X) and y = log(Y ). This has the added advantage that the concentrations X and Y are automatically positive (however Q and F are not necessarily positive), and has become the standard way to solve the problem (see [START_REF] Van Der Lee | CHESS, another speciation and surface complexation computer code[END_REF] for more details).

If we rewrite our chemical system [START_REF] Dzombak | Surface complexation modeling-Hydrous ferric oxide[END_REF] in terms of the new unknowns x and y, we get

       [ log(C) log(S) ] = A [ x y ] + [ log(K C ) log(K S ) ] [ T W ] = A T [ C S ] + [ exp(x) exp (y)
] that can be rewritten as

A T exp ( A [ x y ] + [ log(K C ) log(K S ) ]) + [ exp(x) exp (y)
] -

[ T W ] = 0 (17)
The chemical system [START_REF] Dzombak | Surface complexation modeling-Hydrous ferric oxide[END_REF] can be now seen as searching for a root of the function defined in [START_REF] Kraütle | General multi-species reactive transport problems in porous media: Efficient numerical approaches and existence of global solutions[END_REF]. Numerically, our chemical solver uses the GSL library 1 (see for instance or [START_REF] Galassi | GNU Scientific Library Reference Manual -Revised Second Edition[END_REF]) and KNM library 2 . The chemical solver can use either Newton's method, Newton's method with a finite difference Jacobian, Broyden's method, and the hybridj algorithm (which is a modified version of Powell's Hybrid method, see [START_REF] Powell | A hybrid method for nonlinear equations[END_REF] or the references given for GSL). Due to the way the code is written, user defined algorithm could easily be integrated.

As we have seen above (equation ( 12)), in the context of reactive transport what is of interest is not so much the individual concentrations, but the way they are distributed between their aqueous part Q and their fixed part F . Accordingly, after system [START_REF] Kraütle | General multi-species reactive transport problems in porous media: Efficient numerical approaches and existence of global solutions[END_REF] has been solved, we compute F as in equation ( 12):

F = A T 3 exp(A 3 log(X) + A 4 log(Y )).
To formulate the coupled problem, it will be convenient to condense the solution of the chemical sub-problem by the function Ψ C defined as

F = Ψ C (T ). (18)
It should be emphasized that the function Ψ C is defined implicitly, via the solution of the nonlinear system [START_REF] Kraütle | General multi-species reactive transport problems in porous media: Efficient numerical approaches and existence of global solutions[END_REF]. Also, we have omitted the dependence on W as, in this work, it is assumed constant in time. The "role" of the function Ψ C is thus to take a given total concentration, and to extract its fixed concentration, according to the chemical equilibrium, This will enable us to treat the chemical solver as a black-box.

Transport solver

We want to be able to couple different transport modules. To achieve this goal, we will formulate the result of a transport step in a way analogous to the chemical solver in section 2.1.2. We write the transport problem, after discretization in space as

dZ dt + L(Z) = r,
where r takes into account sources and boundary conditions. We assume that the result of a transport step (after discretization both in time and space) can be put in the form

Z n+1 = Ψ T (Z n , r n , r n+1 ).
For instance, if we use a method of lines type of discretization, with a backward Euler scheme in time, we obtain

Z n+1 = (I + ∆tL) -1 (Z n + ∆tr n+1 ).
The coupling algorithm in section 2.3 will only use the existence of the function Ψ T , not its specific form.

To be specific, we now give the details for the discretization of the transport step used in this work. Rather than the full backward Euler time discretization shown above, we use a splitting technique developed by Siegel et al. [START_REF] Siegel | Solution of the advection-dispersion equation using a combination of discontinuous and mixed finite elements[END_REF] (see also [START_REF] Mazzia | A time-splitting technique for the advection-dispersion equation in groundwater[END_REF]), that enables us to use different spatial and temporal approximations for different transport processes (advection and diffusion). As explained in these references, it is important to be able to use different methods for the discretization of the different phenomena involved in the transport operator. Explicit methods are better suited for dealing with advection type phenomena (implicit methods add more numerical diffusion than explicit ones), whereas diffusion phenomena require an implicit method to overcome the stringent stability restriction of explicit ones. The method is based on operator splitting: each time step has an (explicit) advection half-step, followed by an (implicit) diffusion half-step. We first solve, by an explicit method

Z n+1/2 -Z n ∆t + ∇.(ωuZ n ) = 0 in Ω (19)
Then we solve

Z n+1 -Z n+1./2 ∆t -∇.(D(u)∇Z n+1 ) = 0, in Ω. ( 20
)
Let us detail each half-step: Advection step: We discretize equation ( 19) by an upwind, cell-centered finite volume method. We denote by ZK the average of the concentration over grid cell K. We integrate the equation over a grid cell K, and use Green's formula to obtain

|K| Zn+1/2 K -Zn K ∆t + ∫ ∂K ω u.n Z in K dγ = 0,
where Z in K is the upwind concentration across the boundary of the cell K. Across one face F = K ∪ K′ , the upwind concentration is defined by

Z in K =      Z n K if u.n F > 0 Z n K ′ if u, n ≤ 0 and face F is not a Dirichlet face, Z dF if face F is a Dirichlet face.
Diffusion step: We use the same mixed finite element method as for the flow problem in section 2.1.1.

Indeed, the diffusion half-step can be written in the form

     Zn+1 K - Zn+1/2 K ∆t + ∇.j = 0 j = -D(u)∇ Zn+1 K .
in Ω, which has the same form as [START_REF] Galassi | GNU Scientific Library Reference Manual -Revised Second Edition[END_REF] above, up to the zeroth order term Zn+1/2 K , and this can easily be included in a mixed formulation.

The explicit advection step has a CFL restriction that takes the form

∆t ≤ min K ω|K| max {K,u.n<0 on ∂K }(U K , 0) . ( 21
)
A further refinement is to take several advection steps per diffusion step, by choosing a large diffusion time step, and restricting the advection time step to obey the CFL limit given in [START_REF] Marchand | Analyse de sensibilité déterministe pour la simulation numérique du transfert de contaminants[END_REF].

When we apply this solution method to the coupled problem [START_REF] Feike | Analytical solution of the one dimensional advection equation and two-and three-dimensional dispersion equation[END_REF], the unknown function is Q n+1 , and the right hand side r depends on boundary conditions (known) and both F n (known) and F n+1 (unknown in principle, but given in the context of a fixed point algorithm). It will be convenient to abuse notation, and to denote the solution of the transport step for equation [START_REF] Feike | Analytical solution of the one dimensional advection equation and two-and three-dimensional dispersion equation[END_REF] by

Q n+1 = Ψ T (F n+1 ) (22)
Of course the function Ψ T used in [START_REF] Mazzia | A time-splitting technique for the advection-dispersion equation in groundwater[END_REF] may depend on n.

The code for implementing both the flow solver described in section 2.1.1 and the transport solver described here is based on a the LiveV finite element library 3 .

Solution of the coupled system by a fixed point method

We recall the definition of Ψ C (18) and Ψ T [START_REF] Mazzia | A time-splitting technique for the advection-dispersion equation in groundwater[END_REF], and add the relation

T = Q + F to obtain the discrete coupled system      Q n+1 = Ψ T (F n+1 ) T n+1 = Q n+1 + F n+1 F n+1 = Ψ C (T n+1 ) (23)
This system can be solved by various methods. The classical way is still based on a fixed-point method [START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF] or [START_REF] Mosé | Efficiency of operator splitting procedures for solving reactive transport equation[END_REF]. Recently, more global approaches have been used with good success, based on variants of Newton's method. Several papers implementing global methods are [START_REF] Amir | Newton -Krylov methods for coupling transport with chemistry in heterogeneous porous media[END_REF][START_REF] De Dieuleveult | A numerical model for coupling chemistry and transport[END_REF][START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physics-based preconditioning to biogeochemical transport[END_REF][START_REF] Kraütle | A new numerical reduction scheme for fully coupled multicomponent transport-reaction problems in porous media[END_REF][START_REF] Kraütle | A reduction scheme for coupled multicomponent transport-reaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions[END_REF]. The paper [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF] explains why the so called fixed-point method should more properly be viewed as a Block Gauss-Seidel method.

In this work, we will present results for a fixed-point method, as this has the same ingredient as the more tightly coupled methods, and provides a base case for validating the latter methods. The main advantage of the fixed point method is that the sub-problems Ψ T and Ψ C are decoupled, and this lets us reuse the solvers for the sub-problems that we detailed earlier, namely the chemistry solver described in section 2.1.2, and the transport solver described in section 2.1.3. We note that a Newton-Krylov method for solving the coupled problem ( 23) is currently being implemented. In this method, the nonlinear system ( 23) is solved by Newton's method, and the linear system arising at each Newton iteration is solved by an iterative method, usually GMRES. The main advantage of this method is that it does not require storing the Jacobian matrix of the coupled system, which could become quite large for realistic 3D simulations, with a large number of chemical species. We do not report results for this method, noting that most of the modules developed for the fixed-point method can be reused for Newton's method, but in addition modules for the derivatives are needed. The Newton-Krylov method has been used in the same context by Hammond et al. [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physics-based preconditioning to biogeochemical transport[END_REF], and with the formulation described here, but with only a 1D geometry by Amir and Kern [START_REF] Amir | Newton -Krylov methods for coupling transport with chemistry in heterogeneous porous media[END_REF]. It is this work that is currently being extended to the setup described in this paper.

Implementation of the fixed point method

Suppose that Q n , F n and T n are known. We start by initializing F n+1,0 = F n , and we loop over the following algorithm :

(1) compute Q n+1,k+1 = Ψ T (F n+1,k ) by using the transport solver, once for each aqueous species, (

) compute T n+1,k+1 = Q n+1,k+1 + F n+1,k , ( 2
) compute F n+1,k+1 = Ψ C (T n+1,k+1 3
) by using the chemical solver, once for each degree of freedom of the mesh (4) start again from (1) until the norm ∥F n+1,k+1 -F n+1,k ∥ is small enough (in our test case, small enough means 10 -6 in absolute norm).

Step (2) may look troublesome, as we add aqueous and fixed concentrations from different iterations. Indeed, this is probably the main reason why the so called "Standard Non Iterative Algorithm" (just one iteration of the above algorithm) fails. However, since we iterate to convergence, this discrepancy should not matter as long as the tolerance is tight enough. To set up the initial concentrations, we consider that T 0 and , W are piecewise constants. By using the chemical solver, we compute F 0 = Ψ C (T 0 ) and then Q 0 = [T 0 , W ] -F 0 for each subdomain. 3 LifeV is the joint collaboration between three institutions: École Polytechnique Fédérale de Lausanne (CMCS) in Switzerland, Politecnico di Milano (MOX) in Italy and INRIA (REO) in France, see http://www.lifev.org

Numerical results

Validation of the chemical solver

We report on the performance of our chemical solver on several examples taken from the literature.

Ion exchange example

The first example is used as a validation test case. It comes from the documentation of the PhreeqC [START_REF] Parkhurst | User's guide to PHREEQC (version 2)-A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations[END_REF] software, a standard geochemical code, where it is known as "Example 11". The system is an ion exchange, where a solution of sodium and potassium in equilibrium with an exchanger is injected with calcium-chloride. The reactions involved are as follows

Na + + X -⇌ NaX K + + X -⇌ KX 1 2 Ca 2+ + X -⇌ 1 2 CaX 2
where NaX, KX and CaX 2 are complexes formed on the surface of the rock matrix, and X represents an exchange site with charge -1 (see [START_REF] Dzombak | Surface complexation modeling-Hydrous ferric oxide[END_REF] for details on modeling surface reactions). The corresponding Morel tableau is given by For this chemical system, one can see that there are no secondary aqueous species, i.e. we have Q = X (X

K + Cl -Ca 2+ Na + X - K CaX 2 0 0 1 0 2 2.86 10 3 NaX 0 0 0 1 1 1 KX 1 0 0 0 1 5.01 (24)
Table 1. Morel tableau for ion exchange example here is the fictitious concentration of the exchange site, not it's name).

The quantity that may be determined experimentally is Q and not T . If we start with X and Y (here, we start with X = [2 10 -4 ; 0; 0; 10 -3 ] and Y = [0]), using respectively the Mass Action Law (first line of 12) and the Conservation Law (second line of ( 12)), we may compute T and W (here, T = [7.51 10 -4 ; 0; 0; 1.55 10 -3 ] and W = [1.10 10 -3 ]). We can then recompute

[X app , Y app ] = [T, W ] -Ψ C ([T, W ]) (since Q = X)
and match up X app with X (and Y app with Y ). We compare different algorithms, as implemented in the GSL library. For each algorithm, we start with the initial guess [x, y] = max(log 10 [T, W ], -10) (of course we start with x i = -10 if associated T i ≤ 0). We run each algorithm with a relative tolerance of 10 -8 and an absolute tolerance of 10 -5 , and limit the number of iterations to 2000 (considering the algorithm failed if this limit number is reached).

We have compared the classical Newton's method, without any line search, with two variants of Powell's hybrid method [START_REF] Powell | A hybrid method for nonlinear equations[END_REF], the difference being in the scaling used for the unknowns, and the fact that the hybrids method uses a finite difference Jacobian. The results are shown in table 2.

Algorithm used In table 2, #f and #j are respectively the number of function evaluations and the number of Jacobian evaluation. The Newton algorithm needs the least number of function evaluations, while the hybridj algorithm uses more function evaluations, but only 3 Jacobian evaluations. The finite difference algorithm is not competitive with the methods based on an exact Jacobian. Finally, we note that the hybridj algorithm is more robust, as it converges for several examples when the Newton algorithm fails -see the second test case described in sub-section 3.1.2. In fact, for X Cl -, X Ca 2+ and Y X -, one can never have exactly zero (because we solve for the logarithms), but the result obtained is so much smaller than the relative tolerance that we can consider it is effectively zero.

X K + X Cl -X Ca 2+ X Na + Y X -Iters #f #j Newton 2.

MoMaS Benchmark (a synthetic example)

The first test case involved only ion exchange. We test now the chemical solver on a more complex system, the "easy test case" of the MoMaS benchmark [START_REF] Carrayrou | presentation of the MoMaS reactive transport benchmark[END_REF], and see also [START_REF] Momas | Exercices de qualification numérique de codes[END_REF]. Indeed, one should note that "easy" refers to the chemical phenomena taken into account (here no minerals, and no kinetics), as this test case is by no means easy numerically. We recall here the Morel tableau for the chemical system. As can be seen on the tableau, the equilibrium constants span more than 45 orders of magnitude, and some of the stoichiometric coefficients are quite large, making the system very nonlinear. We take for total concentrations T = [0, -2, 0, 2, 1] and either

X 1 X 2 X 3 X 4 Y K C 1 0 -1 0 0 0 10 -12 C 2 0 1 1 0 0 1 C 3 0 -1 0 1 0 1 C 4 0 -4 1 3 0 0.1 C 5 0 4 3 1 0 10 35 S 1 0 3 1 0 1 10 6 S 2 0 -3 0 1 2 0.1 (25)
Table 3. Morel tableau for MoMaS example W = 1 (test case A) or W = 10 (test case B). We compute the concentrations by using the same algorithms as in sub-section 3.1.1. The method was validated against a solution obtained by using a Gröbner basis method for solving the polynomial system with a symbolic algebra system. The results are shown in table 4. One can see on this table that the "pure" Newton's method fails, whereas both the hybridj algorithm and the hybrids algorithm converge.

Algorithm used Iters

(A) #f (A) #j (A) Iters (B) #f (B) #j (B)

A difficult aqueous system

To see the effects of the initial point choice, we test with another chemical system from [START_REF] Brassard | A feasible set for chemical speciation problems[END_REF], or [START_REF] Carrayrou | New Efficient Algorithm for Solving Thermodynamic Chemistry[END_REF]. We recall the Morel tableau of this system (which contains only aqueous species) on table 5.

We take for total concentrations [T ] = [0, 10 -3 , 10 -3 ], and for both newton and hybridj algorithm we take as initial guess with [x] = [-10, α, β] (in logarithmic form), with α and β varying from -12 to -2, and we limit the maximum number of iterations to 200. On figures 2 and 3, we plot number of iterations needed for convergence with respect to α (horizontal axis) and β (vertical axis). One can see by comparing figures 2 and 3 that for many choices of the initial guess, Newton's method is more efficient that the hybridj algorithm. However, we can see a zone (dark red on figure 3) for which the Newton algorithm fails -however the hybridj algorithm still converges.

H + Al 3+ H 3 L log 10 K OH - -1 0 0 -14 H 2 L - -1 0 1 -4.
These tests tend to support the (well known) conclusion that the pure form of Newton's method is not robust, but that some form of globalization (the hybridj method uses a trust region globalization) makes it much less susceptible to failure when starting far from the solution. This is an especially important issue when the chemical solver is used as a module in a coupled code, as it will be called repeatedly, with little possibility of restarting if it fails. up of two media : medium A (diagonally dashed) is less permeable and less dispersive than medium B (the barrier, horizontally dashed). The characteristics for both media are given in table 6 (taken from the MoMaS benchmark description [START_REF] Carrayrou | presentation of the MoMaS reactive transport benchmark[END_REF][START_REF] Momas | Exercices de qualification numérique de codes[END_REF]).

Medium A Medium B Permeability K (LT -1 ) 10 -2 10 -5 Porosity ω (-) 0.25 0.5 Dispersivity α L (L) 10 -2 6 10 -2 Dispersivity α T (L) 10 -3 6 10 -3 Table 6. Characteristics for the MoMaS benchmark flow and transport example Note that the contrast in permeabilities for the two media is 1000, so this model presents a fairly stringent test for the ability of the code to handle heterogeneous media. The boundary conditions for the flow computations are as follows:

• imposed flux on the boundaries labeled "In" ωu 1 = ωu 2 = 2.25 10 -2 LT -1 ,

• imposed head on the boundary labeled "Out" H = 1L,

• no flow on the remaining boundary and the boundary conditions for transport are

• given concentration on the inflow "In" boundaries Q = 1

• zero diffusive flux on the outflow "Out" boundary,

• zero total flux on the remaining boundary -this condition may be decomposed into two terms :

no convective flux because we have homogeneous flow velocity condition, no diffusive flux, giving an homogeneous Neumann condition, and we take zero initial concentration at time t = 0. We computed the flow with using the LiFeV environment, and saved it as input for the transport solver.

Coupled system

We now consider the full coupled problem, with the chemical system given by [START_REF] Montarnal | Coupling methodology within the software platform alliances[END_REF], subject now to both transport and chemistry. We consider the simple geometry shown on figure 7: We assume that the the flow gives rise to a uniform horizontal velocity field given by u = (1, 0) T . We thus do not need to compute the flow in this case.

With the geometry defined above, we consider two different configurations, depending on the way the CaCl 2 is injected: in this first case, the injection occurs along the full left side of the domain, leading to a solution depending only on x (1D configuration). In the second case, the injection occurs only on part of the boundary, giving rise to a truly 2 dimensional solution (2D configuration).

1D configuration

In this case, the concentration is prescribed along the full left boundary of the computational domain. This corresponds to the ion exchange example from the PhreeqC documentation a;ready mentioned in section 3.1.1. This reference solution will enable us to validate the coupled solution.

For the total aqueous concentrations Q j , we consider :

• Q = [ 0; 1.2 10 -3 ; 6 10 -4 ; 0] on Γ r ∪ Γ g , • ∇Q • n = 0 on Γ b ∪ Γ k .
In other words, we start from a porous media which contains K + and Na + , and we inject on the left side of this media CaCl 2 , letting the right side free. The idea is, since the Darcy velocity u is parallel to the x-axis, to get concentrations T , Q and F that only depend on x (assuming that no there is no transverse diffusion, the transport phenomena and boundary conditions are invariant with respect to y). • Q Cl -(blue line) : since chlorine does not react in our chemical system, it is only subject to transport. The obtained curve coincides with the well-known solution for a perfect tracer. • Q Ca 2+ (black line) : if calcium were not affected by chemical reactions, at each time we would have Q Ca 2+ (t) = 0.5Q Cl -(t) on the right boundary (since the initial data and boundary conditions satisfy this relation). However, calcium does react, and the reaction is quasi-total (since the equilibrium constant K associated to the reaction Ca 2+ + 2X -= CaX 2 is much larger than 1). This explains why calcium reaches the right boundary so late. • Q K + (green line) and Q Na + (red line) : since Ca 2+ reacts with X -, we see a disappearance of KX and NaX. This effect would increase Q K + and Q Na + . However, we have also transport-diffusion effects, and since we inject no sodium nor potassium, these effects will decrease Q K + and Q Na + . So there are two phenomena in competition : transport-diffusion on the one hand, and chemistry on the other hand. For sodium, transport-diffusion is always more important than chemistry. For potassium,the chemistry effects dominate at first, then the transport takes over.

2D configuration

In this second case, we only inject the CaCl 2 solution on part of the left boundary. For the total aqueous boundary conditions, we consider the following conditions :

• Q = [ 0; 1.2 10 -3 ; 6 10 -4 ; 0] on Γ r ,

• ∇Q • n = 0 on ∂Ω \ Γ r . The object of this test is to show truly 2D effects, and not only a 1D solution computed on a 2D grid. We show on figure 10 the aqueous concentration of Cl -(top-left corner), Ca 2+ (top-right corner), K + (bottom-left corner) and Na + (bottom-right corner) at time t = 35. We find results that can be explained according to the theory. In particular, chlorine is still a perfect tracer, and an analytical solution has been computed by Feike and Dane [START_REF] Feike | Analytical solution of the one dimensional advection equation and two-and three-dimensional dispersion equation[END_REF], and our results show a good agreement with this solution. The other concentrations can be explained from the chemical point of view as in the previous section.

Figure 1 .

 1 Figure 1. The Raviart-Thomas space on a rectangle

Figure 2 .Figure 3 .

 23 Figure 2. Convergence of hybridj algorithm

Figure 4 .

 4 Figure 4. Geometry of the domain Ω for the transport-diffusion part

Figure 5 .

 5 Figure 5. Plot of Darcy velocity and pressure Figure 6. Plot of Q Cl - We superimpose on figure 5 the level lines for the pressure on top of the Darcy velocity, and we plot on figure 6 the level lines of the tracer concentration at time t = 25. Qualitatively, we can see that concentration seen on figure 6 follows the Darcy velocity shown on figure 5.

Figure 7 .

 7 Figure 7. Geometry of the domain Ω for the full test

Figure 8 .Figure 8

 88 Figure 8. Plot of different concentrations Q j for time t=20

Figure 9 .

 9 Figure 9. Plot of different aqueous concentrations as functions of time

Figure 10 .

 10 Figure 10. Instants snapshots of the aqueous concentration of Cl -(top-left corner), Ca 2+ (top-right corner), K + (bottom-left corner) and Na + (bottom-right corner)

Table 2 .

 2 Comparison of algorithms for the chemical problem -first test case

		00 10 -4	0	0	1.00 10 -3	0	20	20 20
	hybridj	2.00 10 -4	0	0	1.00 10 -3	0	35	35 3
	hybrids	2.00 10 -4	0	0	1.00 10 -3	0	52	57 0

Table 4 .

 4 Comparison of algorithms for the chemical problem -MoMaS benchmark test case

	Newton	fail			fail
	hybridj	114	114	7	109	109	4
	hybrids	77	97	0	102	137	0

Table 5 .

 5 Morel tableau for the gallic acid example

	15

Acknowledgments

Part of the work reported in this paper took place during the Cemracs summer school in August 2008. The authors thank the organizers of the Cemracs for making possible a very productive (and pleasant) session.

* This work was supported by ANR SHPCO2 project, and partly supported by Groupement MoMaS (PACEN/CNRS, ANDRA, BRGM, CEA, EdF, IRSN)


Lambiénou Yé
email: ylambienou@yahoo.fr

Luc Abbadie

Gérard Bardoux

Jean-Christophe Lata

Hassan Bismarck Nacro

Dominique Masse

Henri De Parseval

Sébastien Barot

Contrasting impacts of grass species on nitrogen cycling in a grazed Sudanian savanna

Keywords: 15 N, nitrogen cycling, savanna, perennial grass, annual grass, herbivore, Burkina Faso

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

Introduction

Plants influence nutrient cycling and their long term availability [START_REF] Chapman | Plants actively control nitrogen cycling: uncorking the microbial bottleneck[END_REF][START_REF] Hobbie | Effects of plant species on nutrient cycling[END_REF], and this feedbacks on their own growth and primary productivity. This is achieved through many mechanisms [START_REF] Hobbie | Effects of plant species on nutrient cycling[END_REF][START_REF] Knops | Mechanisms of plant species impacts on ecosystem nitrogen cycling[END_REF]. This allows them to cope with abiotic constraints of ecosystems such as the lack of mineral nutrients and contrasting soil types and with biotic constraints such as herbivory or competition between plants [START_REF] Craine | Resource strategies of wild plants[END_REF][START_REF] Grime | Plant strategies, vegetation processes, and ecosystem properties[END_REF]. Among all possible strategies, the distinction between perennial and annual plants is particularly relevant in terms of nitrogen cycle and competitive ability [START_REF] Grime | Evidence for the existence of three primary strategies in plants and its relevance to ecological and evolutionary theory[END_REF][START_REF] Grime | Plant strategies, vegetation processes, and ecosystem properties[END_REF]. Annuals are viewed as opportunistic. They need to acquire their resources as quickly as possible and benefit from the immediately available mineral resources.

Perennials are more likely to influence local nutrient cycles and to benefit from such modifications because their root system, their root exudates and their litter (roots and leaves) interact for a longer period with the same patch of soil and soil micro-organisms living in this patch (Vinton and Burke, 1995). Perennials grow slower than annuals but may benefit from sparing soil mineral nutrient resources, from increasing the local availability of these resources and from decreasing nutrient losses. It has for example been shown in tropical savannas that many perennial bunch grasses are able to inhibit nitrification through the release of particular molecules from their root systems [START_REF] Lata | Grass populations control nitrification in savanna soils[END_REF]Subbarao et al., 2007a;Subbarao et al., 2007b), which increases their own biomass [START_REF] Boudsocq | Modelling approach to analyze the effects of nitrification inhibition on primary production[END_REF]. As these grasses live for several tens of years, they build dense root systems where the close proximity between dead and live roots is likely to induce a very efficient nutrient recycling and reduce nitrogen losses [START_REF] Abbadie | Nitrogen dynamics in the soil-plant system[END_REF][START_REF] Abbadie | Independance of savanna grasses from soil organic matter for their nitrogen supply[END_REF]. In this context, our primary goal was to compare, from the point of view of nitrogen cycling, two annual and two perennial bunch grasses growing in the same West African savanna.

To meet this goal, we measured the natural abundance in 15 N of the biomass of these four grasses and of the soil penetrated by their root system. Indeed, 15 N abundance is considered as a good indicator of ecosystem functioning [START_REF] Dijkstra | N-15 enrichment as an integrator of the effects of C and N on microbial metabolism and ecosystem function[END_REF][START_REF] Robinson | [END_REF]Staddon, 2004;Templer et al., 2007). First, 15 N abundance is used to determine the origin of soil organic matter [START_REF] Boutton |  13 C values of soil organic carbon and their use in documenting vegetation change in a subtropical savanna ecosystem[END_REF][START_REF] Danso | Misconceptions and practical problems in the use of 15 N soil enrichment techniques for estimating N 2 fixation[END_REF]. Similarly, the 15 N of biomass gives hints on the origin of the nitrogen. For example legumes, due to symbiotic nitrogen fixation, tend to have lower concentrations in 15 N than other plants and thus to decrease soil 15 N concentration through their litter. Second, the isotopic fractionation due to the lower reactivity of 15 N relatively to 14 N allows tracing nitrogen fluxes in the soil. For example, rapid nitrogen cycling and nitrification are thought to increase 15 N abundance [START_REF] Dijkstra | N-15 enrichment as an integrator of the effects of C and N on microbial metabolism and ecosystem function[END_REF], while tighter nitrogen cycling (e.g. less leaching and denitrification) should lower 15 N abundances. The difference between plant and soil 15 N abundances gives hints on plant preference for ammonium vs. nitrate [START_REF] Kahmen | Foliar  15 N values characterize soil N cycling and reflect nitrate or ammonium preference of plants along a temperate grassland gradient[END_REF]. Similarly, the discrepancy between the 15 N signatures of the biomass of an African perennial bunch grass and bulk soil (lower 15 N in the biomass than in the soil organic matter) [START_REF] Abbadie | Independance of savanna grasses from soil organic matter for their nitrogen supply[END_REF] suggested that this grass fosters an efficient nitrogen cycling through the recycling of the nitrogen contained in the roots of its dense and localized root system. Taken together, we predict that our two perennial and two annual grasses should have different 15 N signatures, which should in turn impact soil signatures. For the same reasons the universal positive correlation between leaf 15 N signature and leaf N concentration (Craine et al., 2009b) should also be different between annual and perennial biomasses.

Besides plant influence on nitrogen cycling, other factors such as herbivory and soil properties are likely to impact these cycles, which should in turn modify 15 N isotopic signatures. Cattle or wild herbivores indeed quicken carbon and nutrient cycling by increasing the turnover of the biomass. Herbivore impact on nutrient cycling has been extensively studied. On the one hand, an important issue is to determine whether herbivores could increase the efficiency of nitrogen cycling (decrease nitrogen losses), which could increase primary production on the long term [START_REF] De Mazancourt | Grazing optimization and nutrient cycling: when do herbivores enhance plant production?[END_REF]. On the other hand, they have been shown to impact many fluxes of mineral nitrogen but the direction of these effects is likely to be case dependent [START_REF] Mcnaughton | Promotion of the cycling of diet-enhencing nutrients by African grazers[END_REF]Wardle et al., 2001). For example, cattle have been shown to increase nitrification, denitrification and free nitrogen fixation in a temperate grassland [START_REF] Patra | Effects of management regime and plant species on the enzyme activity and genetic structure of N-fixing, denitrifyning and nitrifying bacterial communities in grassland soils[END_REF]. In turn, the modification of these fluxes impacts 15 N signatures and may either increase [START_REF] Frank | Effects of native grazers on grassland N cycling in Yellowstone National Park[END_REF] or decrease [START_REF] Frank | Ungulate stimulation of nitrogen cycling and retention in Yellowstone Park grasslands[END_REF] soil 15 N abundances. Similarly, soil properties such as soil depth, soil compaction, soil texture or soil organic matter content are likely to impact biomass and soil 15 N signatures [START_REF] Abadín | Relationships among main soil properties and three N availability indices[END_REF] through their effects on soil capacity to retain mineral nutrients, or through their effects on microbial biomass and activities.

Taking into account these interactive impacts of plant species and soil type/herbivores on nitrogen cycling and detecting these interactions using 15 N signatures is interesting in two contrasted ways. On the one hand, such complex interactions could strongly increase the variability in isotopic signatures, which could blur potential effects of targeted treatments (here the plant species). Documenting these interactions and their effects on isotopic signatures is thus methodologically useful to help designing and interpreting other experiments or field samplings. On the other hand, isotopic signatures might allow analyzing the impact of herbivores and different soil types on nitrogen cycling and might help disentangling the interactions between herbivores or soil types and plant species.

For these reasons, we compared two annual and two perennial grasses in a West African savanna where both soil depth and cattle pressure are varying. Taken together, we make the following predictions. If perennials are able to foster a more efficient nitrogen cycling, (1) there should be lower 15 N abundances in their biomass and soil and (2) the relation between 15 N signature and N concentration should be different between annual and perennial leaf biomass. If soil depth and cattle pressure significantly modify nitrogen fluxes,

(3) soil depth and cattle pressure should influence 15 N signatures, (4) plant species and soil depth or cattle pressure should affect 15 N signatures in an interactive way, i.e. the impact of each species should depend on soil depth and cattle pressure.

Material and methods

Study site

The study site is the third management unit of the protected forest of Dindéresso (FCD) The study site has been divided in four blocs of approximately 1.5 ha according to the dominance of perennial and annual grasses and the frequency of grazing: blocs 1 and 2 are dominated by annual grasses and are more grazed by cattle during the rainy season while blocs 3 and 4 are dominated by perennial grasses and are less grazed during the rainy season.

Blocs 1 and 2 are indeed next to the main road, so that they are more easily reached by cattle and shepherds, and blocs 3 and 4 are supposed to host an abundant population of tsetse flies during the rainy season so that they are avoided to limit the risk of cattle infection by trypanosomiasis. Blocs 1 and 2 have shallower soils (at most 55 cm deep, they are indurated ferruginous leached soils). Blocs 3 and 4 have deeper soils (at least 105 cm deep, the indurated layer is deeper). Blocs 1 and 2 are contiguous and so are blocs 3 and 4. Blocs 1-2 and 3-4 are separated by a distance of about 2000 m.

Soil and biomass sampling

The study focusses on four of the dominant grass species: Andropogon pseudapricus, Loudetia togoensis (annuals), Andropogon ascinodis and Andropogon gayanus (perennial bunch grasses). The four grasses are cespitosus but the tussocks of the two perennials are much larger than the ones of the two annuals. The four species have C 4 photosynthesis [START_REF] Breman | Manuel sur les pâturages des pays sahéliens[END_REF]. As a control, bare soils were also sampled. In each of the four blocs, 2 or 3 replicates of each plant species or bare soil areas were randomly selected: an individual for perennial bunch grasses, a small patch covered of individuals of the same species for annuals and a small patch deprived of plants for bare soils. In total, 10 replicates were sampled for each species and bare soils and they were spread over the 4 blocs. Taken together, 50 sampling units (plant species or bare soil) were investigated. In each sampling unit a soil sample was collected in June 2011 from the 0-10 cm layer using a core auger (5 cm diameter). This sample was taken in the middle of tussocks for perennial grasses. Soil samples were then air dried in the shade and sieved at 2 mm. Leaves and roots were collected from each annual or perennial grass sampling unit. Roots were collected using a core auger (5 cm diameter) and washed with water. Leaves and roots were oven-dried at 70°C for 72 hours. In total, 40 samples of roots, leaves and soil were analyzed as well as 50 soil samples (40 samples for annual and perennial grasses and 10 for bare soils).

Soil and biomass analyses

Soil and biomass samples were thinly ground and their content in N, C, 15 N measured by EA-IRMS (Carlo-Erba NA-1500 NC Elemental Analyser on line with a Fisons Optima Isotope Ratio Mass Spectrometer). As usually done (Wang et al., 2010), the contents in 15 N was expressed as relative differences in the ratios 15 N/ 14 N between samples and international standards: δ 15 N (‰) = (( 15 N/ 14 N) sample /( 15 N/ 14 N) standard ) -1) 1000,

The international standard for N is the atmospheric N 2 .

Statistics

All variables were analysed using ANOVAs testing for the effect of the bloc, the grass species and the interaction between the two (Tables 1 and2). When the interaction was not significant it was removed from the model. When only simple effects were significant the direction of their effects was determined using the estimated model parameters and post-hoc Tukey tests. When the interaction was significant the variable was reanalyzed separately for each grass species (or bare soil). Graphs (Figure 1-2) only display results for factors and factor combinations that are significant. The relation between biomass δ 15 N and N concentration was studied using ANCOVAs testing for the effect of N concentration, the effect a categorical variable (root vs. leaves or the grass species) and the interaction between the two. All these analyses were achieved using R software (R development core team, 2010).

Results

Soil δ 15 N is significantly lower (Tab. 1 and Fig. 1) for the two perennial species (respectively 3.98 and 3.88 ‰ for A. ascinodis and A. gayanus, see Table 3 for all values) than for the two annual species and the bare soil (respectively 4.38, 4.55 and 4.56 ‰ for A. pseudapricus and L. togoensis and bare soil), but is not affected by the bloc. This corresponds to a decrease in δ 15 N of about 0.57 ‰.

Soil N concentration is affected by the interaction between the grass species and the bloc: under A. ascinodis N concentration is higher in blocs 1 and 3 than in blocs 2 and 4, under A. pseudapricus it is higher in blocs 3 and 4 than in blocs 1 and 2, and under L. togoensis it is higher in blocs 1 and 4 than in blocs 2 and 3. Soil C concentration is affected by the interaction between grass species and the bloc (Fig. 2, Tab. 1): under A. pseudapricus it is higher in blocs 3 and 4 than in blocs 1 and 2, under L. togoensis it is higher in blocs 1 and 4 than in blocs 2 and 3. Biomass δ 15 N is overall higher under the two annual than under the two perennial grasses (Fig . 1 and Table 2). Both for roots and leaves there is a significant effect of species and the bloc on biomass δ 15 N but the interaction between these two factors is not significant (Tab. 2). For roots, the δ 15 N is higher (Fig. 2) in the two annuals (respectively 2.60 and 1.26 ‰ for A. pseudapricus and L. togoensis) than in the two perennials (respectively -0.72 and -0.99 ‰ for A. ascinodis and A. gayanus) and in bloc 1 than in blocs 3 and 4 (respectively 1.44, -0.17 and -0.09 ‰ for blocs 1, 3 and 4). This corresponds to an increase in δ 15 N of about 2.8 ‰ between perennials and annuals. For leaves, the δ 15 N is higher (Fig. 2) in A. pseudapricus (2.74) than in the three other species (respectively 0.67, -0.59 and -0.30 ‰ for L. togoensis, A. ascinodis and A. gayanus) and in blocs 1 and 2 than in bloc 4 (respectively 1.74, 1.18 and -0.53 ‰ for blocs 1, 3 and 4). Overall, the results on soil and biomass δ 15 N support our hypothesis that perennial grasses lead to lower δ 15 N than annuals (first hypothesis in the Introduction). However, these results are at odds with our third and fourth hypothesis.

Both for roots and leaves there is a significant effect of species and the bloc on the biomass N concentration but the interaction between these two factors is not significant (Tab.

2). Root N concentration is higher (Fig. 2) in the two annuals (respectively 4.09 and 3.65 mg g -1 for A. pseudapricus and L. togoensis) than in the two perennials (respectively 2.52 and 3.71 mg g -1 for A. ascinodis and A. gayanus) and in blocs 1, 2, 3 (respectively 3.39, 3.72 and 3.13 mg g -1 ) than in bloc 4 (2.77 mg g -1 ) and higher in bloc 2 than in bloc 3. This corresponds to an increase in N root concentration of about 48 % between perennials and annuals. Leaf N concentration, as δ 15 N, is higher (Fig. 2) in A. pseudapricus (10.55 mg g -1 ) than in the three other species (respectively 6.32, 6.20 and 7.41 mg g -1 for L. togoensis, A. ascinodis and A. gayanus) and in blocs 2 and 3 than in bloc 4 (respectively 8.24, 8.26 and 6.41 mg g -1 for blocs 2, 3 and 4).

A first ANCOVA shows that there is an overall positive effect of the biomass N concentration on the δ 15 N and that the slopes of the relation between the two variables is steeper for roots than for leaves (Fig. 3). A second ANCOVA shows that for roots, the relation is the same for all species (δ 15 N= 1.95 N R -5.64, where N R is the root biomass N concentration). A third ANCOVA shows that for leaves the relation varies significantly between species (Fig. 3): the slope is steeper for L. togoensis (δ 15 N= 1.33 N L -7.74, where N L is the leaf biomass N concentration) than for the three other species (δ 15 N= 0.50 N L -3.78).

These results only partially support our second hypothesis (see the Introduction). We also checked that the bloc does not influence the relation between δ 15 N and N concentration in roots and leaves.

Discussion

Our results support, at least partially, two of the hypotheses (the first and the second) we set in our introduction: (1) 15 N signatures suggest that perennial and annual grasses do not have the same impact on nitrogen cycling. (2) The increase in leaf δ 15 N with N concentration is steeper in L. togoensis than in the three other grasses. (3) 15 N signatures do not allow to point at any overall effect (whatever the grass species) of cattle grazing or soil depth on grass functioning or nitrogen cycling. (4) 15 N signatures do not suggest any interaction between grass species and blocs that would affect nitrogen cycling.

Impact of grasses on δ 15 N

The δ 15 N values are higher in the soil of annual grasses and bare soil than in the soil of perennial grasses and, consistently, the same qualitative pattern is found for the δ 15 N of the grass biomasses. Since the source of nitrogen should be overall the same for plants growing next to each other, this pattern must be due to fine scale differences between species in the various nitrogen fluxes leading to isotopic fractionation. The observed higher soil δ 15 N under annual grasses and in bare soil can be interpreted as an indicator of increased rates in nitrogen cycling (mineralization, nitrification and denitrification) that also foster nitrogen losses [START_REF] Nacro | Natural nitrogen-15 abundance and carbohydrate content and composition of organic matter particle-size fractions from a sandy soil[END_REF]Templer et al., 2007;Templer et al., 2008). Indeed, fractionation likely leads to an enrichment of the residual pool in 15 N because 14 N, the lighter of the two stable isotope, is more reactive [START_REF] Mariotti | Experimental determination of nitrogen kinetic isotope fractionation: Some principles; illustration for the denitrification and nitrification processes[END_REF]. Thus, in the case of increased rates of nitrogen cycling, 14 N is more likely to leave the ecosystem through denitrification and leaching, which leads to an increase in soil δ 15 N. This interpretation is fully compatible with a previous study [START_REF] Abbadie | Independance of savanna grasses from soil organic matter for their nitrogen supply[END_REF]) that compared the 15 N signature of soil below perennial bunch grasses (+3.4 ‰) and between these grasses, i.e. a bare soil situation (about +5 ‰). Here, we found qualitatively the same difference (+4.0 ‰ for perennials vs. +4.5 ‰ for annuals and bare soils) that can be interpreted as in [START_REF] Abbadie | Independance of savanna grasses from soil organic matter for their nitrogen supply[END_REF]. There are differences in nitrogen fluxes and fractionation between the soil below perennial grass tufts and the soil between these tufts or below annuals grasses. These differences allow differences in 15 N signatures to build up slowly along the several tens of years of the life of these perennials. Indeed, the soil below grass tufts gets slightly impoverished in 15 N (relatively to the soil between tufts or below annuals), grasses take up this nitrogen so that the biomass of their root also gets impoverished in 15 N, some roots die and release mineral nitrogen impoverished in 15 N, differences in fractionation further decrease the 15 N signature of the nitrogen that can be absorbed by other roots of the same grass tuft, and so on… Other studies have shown that soil and biomass δ 15 N values depend on plant species [START_REF] Kahmen | Foliar  15 N values characterize soil N cycling and reflect nitrate or ammonium preference of plants along a temperate grassland gradient[END_REF][START_REF] Kriszan | Long-term changes of the δ 15 N natural abundance of plants and soil in a temperate grassland[END_REF][START_REF] Nadelhoffer | 15 N natural abundances and N use by tundra plants[END_REF]Templer et al., 2005;Wang et al., 2010). As hypothesized in the introduction, we suggest that perennial grasses are able to influence the nitrogen cycle in a way that increases the availability of nitrogen in their favor. The fact that δ 15 N values of bare soils are identical to the values of annual grasses supports this rationale. Several mechanisms leading to a lower fractionation below perennial bunch grasses than between these grasses and below annuals are probably leading to this general pattern: (1) Many African perennial bunch grasses (e.g. Hyparrhenia sp., Brachiaria sp.) have been shown to inhibit nitrification [START_REF] Lata | Grass populations control nitrification in savanna soils[END_REF]Subbarao et al., 2007a;Subbarao et al., 2007b). Perennial grasses are able to build a dense and perennial root system so that (2) the proximity between dead and living roots reduces the risk of leaching of the mineral nitrogen released by root decomposition [START_REF] Abbadie | Nitrogen dynamics in the soil-plant system[END_REF][START_REF] Abbadie | Independance of savanna grasses from soil organic matter for their nitrogen supply[END_REF] and (3) rates of leaching are likely to be further decreased by the fact that living roots are always present and able to take up available mineral nutrients [START_REF] Joffre | Plant and soil nitrogen dynamics in mediterranean grasslands: a comparison of annual and perennial grasses[END_REF]. Testing fully our hypotheses and interpreting thoroughly the 15 N signatures will off course require (1) documenting precisely nitrogen fluxes such as nitrification and denitrification potentials, (2) measuring the 15 N signatures of all nitrogen sources, (3) taking into account mechanisms that could complicate the interpretation of 15 N signatures such as nitrogen fixation by free bacteria or endophytes [START_REF] Elbeltagy | Endophytic colonization and in planta nitrogen fixation by a Herbaspirillum sp. isolated from wild rice species[END_REF].

Impact of soil and herbivores on δ 15 N

Herbivores [START_REF] Aranibar | Nitrogen isotope composition of soils, C-3 and C-4 plants along land use gradients in southern Africa[END_REF]Craine et al., 2009a;[START_REF] Frank | Effects of native grazers on grassland N cycling in Yellowstone National Park[END_REF][START_REF] Frank | Ungulate stimulation of nitrogen cycling and retention in Yellowstone Park grasslands[END_REF] and soil properties [START_REF] Abadín | Relationships among main soil properties and three N availability indices[END_REF] have often been shown to impact soil δ 15 N, which is an indication of modifications of nitrogen cycling. In particular, herbivores impact nitrogen cycling and aboveground-belowground linkages both directly through changes in plant growth and physiology and alterations of the quality and quantity of mineral and organic resources and, indirectly, through changes in the functional composition of vegetation [START_REF] Bardgett | Herbivore-mediated linkages between aboveground and belowground communities[END_REF]. Here, the species-specific δ 15 N of the soil does not depend on the four blocs while root and leaf δ 15 N tends to be higher in the blocs with the higher cattle pressure (and shallower soil). This suggests (1) that increasing cattle pressure and decreasing soil depth has a relatively weak direct effect on nitrogen cycling and/or that this effect is recent (effects on biomass δ 15 N but no effect on soil δ 15 N), ( 2) that this effect tends to increase the rate of nitrogen cycling and nitrogen losses [START_REF] Dijkstra | N-15 enrichment as an integrator of the effects of C and N on microbial metabolism and ecosystem function[END_REF]Templer et al., 2007).

Besides, our sampling design is based on measurements made at the scale of individual tufts for perennial grasses and small mono-specific patches for annual grasses and we investigated equally the same four grass species in the four blocks. In fact, grazing pressure and soil depth do not change the impact of each grass species on soil δ 15 N but these factors likely strongly impact nitrogen cycling and 15 N signatures, indirectly, through their effects on vegetation. Perennial grasses are dominant in the two blocs with deeper soils and lower grazing pressures, while annual grasses are dominant in the two other blocs. This would be the hint that cattle and shallow soils slow down or even reverse the normal succession, which is one of the main mechanisms through which herbivores impact ecosystem functioning [START_REF] Bardgett | Herbivore-mediated linkages between aboveground and belowground communities[END_REF]. Indeed, annuals normally start to grow in fallows and are progressively replaced by perennials [START_REF] Bilgo | Chemical and microbial properties of semiarid tropical soils of short-term fallows in Burkina Faso, West Africa[END_REF]Somé, 1996). Thus, overall, at the scale of the whole savanna lower soil depth and higher cattle pressure lead to higher δ 15 N values through a switch in the dominant grass species (see the first section of the Discussion). Note that the direct effects of cattle pressure also lead to an increase in δ 15 N, at least in biomasses. This is overall a mark of a faster nitrogen cycling and higher nitrogen losses. This is in line with published results [START_REF] Frank | Effects of native grazers on grassland N cycling in Yellowstone National Park[END_REF][START_REF] Patra | Effects of management regime and plant species on the enzyme activity and genetic structure of N-fixing, denitrifyning and nitrifying bacterial communities in grassland soils[END_REF] emphasizing cases of negative effects of herbivores on soil fertility. This also suggests that our blocs 1 and 2 are overgrazed while a milder grazing pressure (blocs 3 and 4) do not lead to such negative effects both in terms of (1) the switch form perennial to annuals grasses and

(2) the long-term soil fertility.

We also predicted that the interaction between grazing intensity and grass species should impact nutrient cycling because species-specific effects on nitrogen cycling could be impacted by grazing, e.g. through differences in palatability. For example, perennial grasses could control efficiently nitrogen cycling and limit nitrogen losses only when they are not grazed intensively and have high biomasses. The absence of significant interactions between bloc and grass species for soil and biomass δ 15 N further supports the idea that at a small scale the impact of each grass species and bare soil drives the local nitrogen cycle and not grazing or soil depth.

Impacts of grasses, soil type and herbivory on nitrogen concentrations

The more efficient nitrogen cycling fostered by perennial than by annual grasses and suggested by δ 15 N values could have led to higher soil nitrogen concentrations for perennials than for annuals and bare soils. Perennial grasses have indeed been found to increase nitrogen content within their root systems [START_REF] Jackson | Geostatistical patterns of soil heterogeneity around individual perennial plants[END_REF]Vinton and Burke, 1995). Such a pattern was not found, and instead, soil nitrogen content depends on a complex interaction between plant species and bloc. This pattern is difficult to explain in details. It seems that the likely interactions between plant species impact on nutrient cycling and grazing by cattle or soil depth is more influential for soil N concentration than for soil or biomass δ 15 N. The pattern must be linked to the fact that the likely reduced rate of nitrogen losses under perennials does not necessarily lead to the buildup of the nitrogen stock within the soil because perennial grasses also store a high amount of nitrogen in their root and shoot biomass. Such an apparent discrepancy between fluxes and stocks was observed in Lamto savanna in Ivory Coast where subpopulations of perennial grasses controlling or not nitrification lead to 2-fold differences in grass biomass but not in soil N content [START_REF] Lata | Stable coexistence of contrasted nitrification statuses in a wet tropical savanna ecosystem[END_REF]. Our results on soil carbon concentration confirms the idea that in this savanna perennial grasses improve their own access to nitrogen but do not necessarily increase locally soil C or N concentrations.

The pattern of nitrogen concentration within grass biomass is simpler: in leaves the nitrogen concentration is higher in A. pseudapricus biomass than in the three other species, in roots the nitrogen concentration is higher in the two annuals than in the two perennials. An explanation could be that annual grasses tend to have higher nutrient concentrations than perennials because they favor a quick growth [START_REF] Craine | Community traitscape of foliar nitrogen isotopes reveals N availability patterns in a tallgrass prairie[END_REF][START_REF] Garnier | Carbon and nitrogen content of congeneric annual and perennial grass species: relationships with growth[END_REF]. This is supported by our root results but only partially by our leaf results: only one annual (A. pseudapricus) has a higher leaf N concentration than the two perennials. Nitrogen concentration depends on the availability of mineral nitrogen but is also linked to plant ecophysiological traits such as the capacity to take up nitrogen or nitrogen use efficiency but we so far do not have a clear explanation for the higher value observed for A. pseudapricus leaves.

Relation between δ 15 N and N concentration in biomasses

Our results are consistent with the already described world-wide positive correlation between δ 15 N and N concentration in leaf biomasses (Craine et al., 2009b). In comparison to a global data base, our leaf N concentrations are rather low (below 12 mg g -1 in most cases) and our leaf δ 15 N values are intermediate (between -2 and 4 ‰ in most cases). Our leaf δ 15 values are relatively high taking into account our rather low N concentration values (see Fig. 3c in Craine et al., 2009b). These values are compatible with the global increase in δ 15 N with the mean annual temperature (our study site has a high mean temperature, about 28°C) and the global decrease with the mean annual precipitation (our study site has an intermediate mean annual precipitation, about 1000 mm yr -1 ) (Craine et al., 2009b).

We expected different slopes for the relation between leaf N concentration and δ 15 N between annuals and perennials. Somehow, a species for which N concentration increases quicker with δ 15 N is able to allocate more nitrogen to its leaves for a given rate of nitrogen cycling and a given rate of openness of this cycling [START_REF] Dijkstra | N-15 enrichment as an integrator of the effects of C and N on microbial metabolism and ecosystem function[END_REF]Templer et al., 2007). This species would be more efficient in terms of exploitation of nitrogen than species with N concentrations increasing slower with δ 15 N. We did find a steeper relation between N concentration and δ 15 N for the annual species L. togoensis. However, we found the same relation for the other annual and the two perennials. All factors influencing the δ 15 N, e.g. mycorrhization and the type of mycorrhiza (Craine et al., 2009b), could explain this difference between L. togoensis and the three other species. More information would be needed to interpret this result.

Our results on the relation between root N concentration and root δ 15 N are consistent with the comparison of 90 grass species across four regions of the world [START_REF] Craine | Environmental constraints on a global relationship among leaf and root traits of grasses[END_REF]. This relation parallels the relation between leaf N concentration and leaf δ 15 N but roots have lower N concentrations than leaves so that the δ 15 N increases more steeply with N concentration for roots than for leaves. For roots the relation between N concentration and δ 15 N is the same for the four species while, for leaves, L. togoensis displays a particular relation. This suggests that the particular strategy of this species for N management aims at increasing leaf N concentration and photosynthesis and not at improving the acquisition of soil resources.

Conclusion

Our main conclusion remains that perennial grasses decrease soil and biomass 15 N signatures relative to annual grasses and bare soil. This confirms the usefulness of 15 N as an integrative tool to asses nitrogen cycling [START_REF] Kahmen | Foliar  15 N values characterize soil N cycling and reflect nitrate or ammonium preference of plants along a temperate grassland gradient[END_REF]Templer et al., 2007). This suggests that perennials better control nitrogen cycling through a slower cycling and lower rates of nitrogen losses through leaching and denitrification: perennial grasses would have evolved particular mechanisms to improve nitrogen cycling in their favor, which would ultimately feedback on their growth allowing them to accumulate more biomass. These arguments are further supported by the fact that annual grasses can be viewed as r or ruderal species in comparison to perennial grasses [START_REF] Grime | Evidence for the existence of three primary strategies in plants and its relevance to ecological and evolutionary theory[END_REF][START_REF] Grime | Plant strategies, vegetation processes, and ecosystem properties[END_REF]. In annuals, evolution should have selected strategies that favor the immediate use of mineral resources. Perennials, on the contrary, should benefit from strategies allowing them to make the best use of mineral resources on the long term, which should lead to lower rates of nutrient cycling and lower rates of nutrient losses.

The mechanisms evolved by perennials could be involved in the high primary productivity of Guinean humid savannas [START_REF] Boudsocq | Modelling approach to analyze the effects of nitrification inhibition on primary production[END_REF][START_REF] Gignoux | Biomass cycle and primary production[END_REF] and could increase the primary productivity of Sudanian savannas dominated by perennial grasses (for example, here in Burkina Faso). In turn, our results are consistent with the idea that annual grasses are a mark of soil and vegetation degradation, and that annual grasses tend to replace perennials in less fertile soil conditions and after over-grazing [START_REF] Belsky | Effects of livestock grazing on stand dynamics and soils in upland forests of the interior West[END_REF][START_REF] Burke | Plant-soil interactions in temperate grasslands[END_REF][START_REF] César | L'influence de l'exploitation sur la pérennité des pâturages de savane II. Rôle du système racinaire dans la dégradation du pâturage[END_REF][START_REF] Derner | Does grazing mediate soil carbon and nitrogen accumulation beneath C 4 , perennial grasses along an environmental gradient?[END_REF]Rossignol et al., 2006). Here, because the higher cattle pressure occurs in the area of shallower soils it was a priori not possible to disentangle their possible effects, which should be possible with an experimental approach excluding cattle from some patches. Anyway, our results suggest that soil or biomass 15 N signatures could be used as a diagnostic tool for the degradation of pastures and the fertility of their soils.

28

Figure 1: Boxplots for N content and δ 15 N of the soil, the roots and the leaves as a function of blocs and grass species. Medians are displayed together with 25th and 75th quartiles.

Treatments and combinations of treatments that are displayed correspond to significant effects as determined by the ANOVAs (Table 1). Different letters denote significant differences between treatments. When the interaction between bloc and grass species is significant different letters denote differences between blocs within a grass species. 1, 2, 3, 4, blocs; A.

pseuda stands for A. pseudapricus. 2). Different letters denote significant differences between treatments. When the interaction between bloc and grass species is significant different letters denote differences between blocs within a grass species.

1, 2, 3, 4, blocs; A. pseuda stands for A. pseudapricus.

  Wang, L., D'Odorico, P.,Ries, L., Macko, S.A., 2010. Patterns and implications of plant-soil δ 13 C and δ 15 N values in African savanna ecosystems. Quaternary Research 73, 77-83. Wardle, D.A., Barker, G.M., Yeates, G.W., Bonner, K.I., Ghani, A., 2001. Introduced browsing mammals in New Zealand natural forests: aboveground and belowground consequences. Ecological Monographs 74, 587-614.

Figure 2 :

 2 Figure 2: Boxplots for soil C concentration. Medians are displayed together with 25th and

Fig. 3

 3 Fig. 3 Regression between biomass δ 15 N and N concentrations. The regression is significantly

in the West of the town of Bobo-Dioulasso, at the altitude of 390

	dominant shrub/small tree species: Vitellaria paradoxa, Terminalia laxifolia, Detarium
	microcarpa, Parkia biglobosa, Guiera senegalensis, Combretum nigricans, Gardenia
	ternifolia. The herbaceous layer is dominated by grasses. The main annual grasses are
	Andropogon pseudapricus, Loudetia togoensis, Microchloa indica. The main perennial
	grasses are Andropogon gayanus, Andropogon ascinodis, Hyparrhenia subplumosa,
	Schizachyrium sanguineum and they are all bunch grasses. There are some legumes (Cassia
	mimosoides, Indigofera trichopoda, Zornia glochidiata, Tephrosia pedicelata, Tephrosia
	bracteolate), Cyperaceae (Fimbristylis hispidula) and other forbs (Waltheria indica, Pandiaka
	heudelotii, Spermacoce stachydea, Striga hermonthica). Grass aboveground biomass and
	necromass are burnt each year by bushfires.
	m (11°12.494' north,
	4°24.159' west). The climate is South-Sudanian: there is a wet season from May to October
	and a dry season from November to April. The area is located between the 900 and 1250 mm
	isohyets. 1254 mm of rain fell in 2010 but only 831 mm in 2011. The mean annual
	temperature is 28°C. The whole forest lies on sedimentary rocks and our own study area lies
	on Bobo-Dioulasso sandstone. Soils are tropical ferruginous leached soils indurated or not.
	Vegetation consists in a shrub savanna grazed by cattle. It is characterized by the following

Table 1 :

 1 Analyses of variance of soil N concentration and δ 15 N and C concentration. The effect of the bloc, the grass species and the interaction between the two are tested. F values are given together with an indication of significance. The last line indicates the direction of effects through the ordering of treatments. When the interaction is significant this line gives the results of ANOVAs analysing separately the bloc effect for each species. **, P <0.01; ***, P<0.001; ns, no significant effect; df, degree of freedom; 1, 2, 3, 4, blocs; AA, A. ascinodis; AG, A. gayanus; AP, A. pseudapricus; LT, L. togoensis; BS, bare soil.

		Df Soil δ 15 N	Soil N	Soil C
	Bloc	3	1.24	2.14	0.61
	Species	4	4.34**	0.57	1.04
	BlocXSpecies	12 -	2.68*	2.09**
	R 2		0.33	0.58	0.57
	Direction of		AP, LT, BS	In AA: 1, 3	In AA: ns
	effects		> AA, AG	>2, 4	In AG: ns
				In AG: ns	In AP: 3, 4>1,
				In AP: 3, 4	2
				>1, 2	In LT: 1, 4>2,
				In LT: 1, 4	3
				>2, 3	In BS: ns
				In BS: ns

Table 2 :

 2 Analyses of variance of root and leaf N concentration and δ 15 N. The effect of the bloc, the grass species and the interaction between the two were tested. The interaction was never significant and was thus removed from the model. F values are given together with an indication of significance. The last line indicates the direction of effects through the ordering of treatments. When the interaction is significant this line gives the results of ANOVAs analysing separately the bloc effect for each species. **, P <0.01; ***, P<0.001; ns, no significant effect; df, degree of freedom; 1, 2, 3, 4, blocs; AA, A. ascinodis; AG, A. gayanus; AP, A. pseudapricus; LT, L. togoensis; BS, bare soil.

		Df Root δ 15 N	Root N	Leaf δ 15 N	Leaf N
	Bloc	3	3.75*	14.10***	4.70**	3.10*
	Species	3	13.88*	44.92***	8.79***	16.63***
	R 2		0.74	0.88	0.73	0.70
	Direction of		AP, LT	AP, LT	AP >LT, AA,	AP >LT, AA,
	effects		>AA, AG	>AA, AG	AG	AG
			1 > 3, 4	1, 2, 3>4	1, 2> 4	2, 3 >4
				2>3

Table 3 :

 3 Means and standard deviations (in parentheses) for biomass and soil N concentrations (mg g -1 ) and δ 15 N (‰).

		Soil δ 15 N	Soil N	Root δ 15 N Root N	Leaf δ 15 N Leaf N
	A. ascinodis	3.98 (0.55)	2.09 (2.76)	-0.72 (0.86)	2.52 (0.31)	-0.59 (1.32)	6.20 (1.95)
	A. gayanus	3.88 (0.60)	2.03 (2.21)	-0.99 (0.73)	2.71 (0.47)	-0.30 (0.52)	7.41 (1.17)
	A. pseudapricus	4.38 (0.34)	3.01 (2.37)	2.60 (1.78)	4.09 (0.45)	2.74 (2.06)	10.55 (2.33)
	L. togoensis	4.55 (0.43)	2.84 (2.75)	1.26 (1.80)	3.65 (0.41)	0.67 (1.92)	6.32 (0.85)
	Bare soil	4.56 (0.44)	2.94 (2.27)	-	-	-	-

Acknowledgments. Lambiénou Yé has been financed by the AIRD program BSTD for his

PhD. We thank the LEMSAT laboratory in Dakar for the facilities we were given for soil and biomass analyses. We also thank Roger Kissou and Saïbou Nignan for their efficient help during the soil characterization and installation of the plots.

Rossignol, N., Bonis, A., Bouzillé, J.-B., 2006. Consequence of grazing pattern and vegetation structure on the spatial variations of net N mineralisation in a wet grassland. Applied Soil Ecology 31, 62-72. Somé, N.A., 1996. Les systèmes écologiques post-culturaux de la zone soudanienne (Burkina Ecology 89, 3030-3040.

Vinton, M.A., Burke, I.C., 1995. Interactions between individual plant species and soil nutrient status in shortgrass steppe. Ecology 76, 1116Ecology 76, -1133. .