R P Brinkmann

The plasma-sheath transition in low temperature plasmas: On the existence of a collisionally modified Bohm criterion

Keywords: numbers: 52, 40, Kh

The plasma-sheath transition in stationary low temperature plasmas is investigated for arbitrary levels of collisionality. The model under study contains the equations of continuity and motion for a single ion species, Boltzmann's equilibrium for the electrons, and Poisson's equation for the field.

Assuming that the electron Debye length λ D is small compared to the ion gradient length l = n i / ∂n i ∂x , a first order differential equation is established for the ion density n i as function of the transformed spatial coordinate q = n i dx. A characteristic feature of this novel 'sheath equation' is an internal singularity of the saddle point type which separates the depletion-field dominated sheath part of the solution from the ambipolar diffusion controlled plasma. The properties of this singularity allow to define, in nonarbitrary way, a collisionally modified Bohm criterion which recovers Bohm's original expression in the collisionless limit but remains meaningful also when collisions are included.

A comparison is made with the collisionally modified Bohm criteria proposed by Godyak

1

The subdivision of a gas discharge into the plasma bulk and the sheath was introduced by Langmuir in 1928 [START_REF] Langmuir | Proc. Natl. Acad. Sci[END_REF]. Together with Tonks he also coined the term sheath edge for the separation of the two "rather distinct" zones [START_REF] Tonks | [END_REF]. (In the standard 1D geometry, bulk and sheath are intervals while the sheath edge is a point.) The connection of the two zones was first studied by Bohm who, in 1949, formulated his famous "criterion for a stable sheath" [START_REF] Bohm | The characteristics of electrical discharges in magnetic fields[END_REF]:

To ensure nonoscillatory solutions, ions must enter the sheath from infinity with a minimum speed (the 'Bohm speed') v B = T e /m i . (Here, m i denotes the ion mass and T e the temperature of the electrons in energy units, i.e., k B ≡ 1.)

The Bohm criterion |v i | ≥ v B has vexed the plasma community for more than four decades.

One of its mysteries was that it could not only be derived from Bohm's original sheath equations (collisionless ions plus Poisson's equation), but also from a "presheath model" (collisional ions plus quasineutrality); there, however, from the properties of a singularity [START_REF] Persson | [END_REF].

A systematic reconstruction of the Bohm criterion was given in 1991 by Riemann who investigated the plasma-sheath transition under the condition that the Debye length λ D is small compared to the mean free path λ [5]. Asymptotically matching the sheath (scale λ D )

and the presheath (scale λ) by a transition layer (scale λ 4/5 D λ 1/5 ), he found a solution uniformly valid on all scales and showed that the Bohm criterion is marginally fulfilled.

Riemann's analysis of the regime λ D /λ ≪ 1 motivated research into what happens when the scale ratio is large. It proved difficult: Some researchers who studied the plasma-sheath transition in the presence of collisions found that the Bohm criterion must be modified [6][7][8][9][10].

Others objected: Riemann stressed that "there is no reason and no basis to formulate a new modified Bohm criterion accounting for collisions in the sheath" and that the "various heuristic attempts" are "inconsistent and lead to unreasonable results" [11]. Franklin went along: "There is no such thing as a collisionally modified Bohm criterion" [12]. This manuscript aims to take part in that discussion by offering an alternative perspective.

A standard sheath model will be analyzed, similar to the ones investigated in refs. [6][7][8][9][10][11][12].

As in [11] and [12], the Debye length λ D will be treated as small compared to a certain other scale of the dynamics. However, that 'other scale' will not be the ion mean free path λ but the ion density gradient length l = n i /| ∂n i ∂x |. Employing arguments that are not heuristic, it will be shown that the assumption λ D ≪ l implies -for all ratios λ D /λ -the existence of a certain mathematical structure (a removable singularity) with properties that allow to unambiguously define a 'collisionally modified Bohm criterion'.

As stated, the model under study is quite conventional. A one-dimensional Cartesian geometry is assumed, described by an x-axis pointing from the electrode x E = 0 to the bulk.

The field is related to the ion and electron densities via Poisson's equation, with the electrons obeying Boltzmann's relation with an electron temperature T e :

-ǫ 0 ∂ 2 Φ ∂x 2 = ǫ 0 ∂E ∂x = e n i -n 0 exp eΦ T e . (1)
Ionization is neglected so that the ion flux density Ψ i to the electrode is spatially constant.

The equation of motion takes into account acceleration by the electrical field and friction due to charge exchange and elastic collisions with the neutrals:

n i v i = -Ψ i , (2)
m i v i ∂v i ∂x = eE - |v i | λ m i v i . (3)
The mean free path λ depends on the velocity of the ions. In many gases, it scales like |v i | for small v i (Langevin interaction) and is constant for large v i (charge exchange interaction and 'hard sphere' elastic collisions). In argon, literature results are well reproduced by

λ(v i ) = |v i |/n N σ 0 v 0 2 + v 2 i
where n N is the gas density and σ 0 = 10 -18 m 2 and v 0 = 550 m/s. The system is completed by the conditions of neutrality and transport equilibrium for x → ∞ and the specification of a potential or current at the electrode location x E .

For comparison, it is advantageous to first study the model without any simplifications. Fig. 1 shows solutions which cover the range from a collisionless (p = 0.1 Pa) to a collision dominated sheath (p = 100 Pa). The ion flux is Ψ i = 2.6 × 10 18 /m 2 s, the electron temperature is T e = 3 eV. All curves exhibit a smooth transition from electron depletion to quasineutrality (= from sheath to plasma) without any pecularities at the Bohm speed v B or the density n B = Ψ i /v B . This, of course, is not surprising; the Bohm criterion is absent in the unabridged model and appears only within asymptotic regimes [11,12].

The standard asymptotic regime in the investigation of the Bohm criterion is 'λ D ≪ λ'.

T e m i - Ψ 2 i n 2 i 1 n i ∂n i ∂x = 1 λ Ψ 2 i n 2 i (4)
Solutions of the presheath model are also shown in fig. 1, calculated for the same data set as the exact solutions to which they asymptotically converge. All have a singularity at the density n B = Ψ i /v B and the speed v B , regardless of the pressure and the collision model.

However, it is only for small values of λ D /λ that this singularity (the 'Bohm point') has any relation to the unabridged model: With decreasing λ, the Bohm point moves from the bulk side of the sheath (where quasineutrality applies) to the center (where quasineutrality fails).

The ambipolar field ceases to be a valid approximation and the presheath solutions begin to deviate from their exact counterparts. Riemann's interpretation of this situation was that with increasing collisionality the "Bohm criterion loses gradually its significance" [11].

However, could that be an over-interpretation? Strictly stated, what loses significance is only the 'λ D ≪ λ approach' to the Bohm criterion, for instance Riemann's own asymptotic analysis in λ D /λ and the uniform approximation based on it. That does not prove that there is no other possibility to reconstruct the Bohm criterion, i.e., to unambiguously single out a certain point in the gradual plasma-sheath transition as the 'sheath edge'. If such a point could be constructed without resorting to the assumption λ D ≪ λ, the resulting 'Bohm criterion' would not lose its significance at higher collisionality.

How could that be done? Any conceivable 'Bohm criterion' should mark the transition from the ambipolar diffusion dynamics of the bulk to the unipolar motion in the sheath.

In the λ D ≪ λ 'presheath approach' it describes the singularity which occurs when inertia equals the ambipolar field and the prefactor of ∂v i ∂x vanishes. A finite λ D /λ model, in contrast, should not terminate but merge into a solution for the sheath. Necessarily then n e = n i and the depletion field must enter the picture. In other words, a more general expression for the field is needed which, however, must still be formulated in terms of the local ion density and its spatial derivative: Otherwise, it could not be combined with the ion inertia term.

Seemingly, such a generalized but local representation of the electric field is impossible: Physically, the field is local only in the quasineutral zone where the ambipolar formula holds; in the depletion region it is not local but dependent on the spatial integral of the ion density.

However, this problem can be overcome: What is really required is not physical locality but only mathematical locality, i.e., locality in a transformed coordinate q(x). The key idea is to take as the new coordinate the spatial integral of the ion density itself, q(x) = n i (x) dx.

The field in the depletion region can be then expressed as a linear function on q, while the field in the ambipolar zone remains proportional to a derivative of n i .

In two recent publications [13,14] the present author has spelled out the details of such a field construction: First, the equivalent electron edge s is defined to be the point where the integrated electron charge below is just equal to the net charge density above:

s x E n e (x ′ ) dx ′ = ∞ s n i (x ′ ) -n e (x ′ ) dx ′ .
(

) 5
The Boltzmann-Poisson equation is then approximately solved under the assumption that the transition from electron depletion to quasineutrality -which takes place within a few Debye lengths around s -is fast compared to the gradient length of the ion density l = n i / ∂n i ∂x . Technically, separate expansions are established, respectively, for the unipolar zone x ≪ s, the transition zone |x -s| < ∼ λ D , and the ambipolar zone x ≫ s. The expansions are systematically carried out up to the order where the ambipolar field appears in the zone x ≫ s, and then spliced together by the 'asymptotic matching principle'. The last step of the construction is an approximate transformation, again employing λ D ≪ l, of the field expression into the 'charge coordinates' defined as

q(x) = x s n i (x ′ ) dx ′ . (6)
The resulting formula for the field, termed the advanced algebraic approximation (AAA), can be seen as an improvement of Godyak's well-known step model [15]. (It reduces to that model in the limit T e → 0.) The AAA is a mathematically local representation of the field in terms of the charge coordinate q, the ion density n i , and the derivative ∂n i ∂q :

E q, n i , ∂n i ∂q = -Ξ 0 eq √ ǫ 0 T e n i T e n i ǫ 0 -Ξ 1 eq √ ǫ 0 T e n i
T e e ∂n i ∂q .

The uniquely defined functions Ξ 0 and Ξ 1 (see fig. 2) are smooth (in fact, analytical) functions of their argument ξ which is the q-version of the Debye length-weighted distance from the location to the equivalent electron edge:

ξ = eq √ ǫ 0 T e n i ≡ q -0 λ D n i ≈ x -s λ D . (8)
Physically, Ξ 0 and Ξ 1 act as 'switches' between the depletion and the quasineutral regime.

The limit q ≪ 0 gives Ξ 0 = -ξ and Ξ 1 = 0 and describes the depletion field; the limit q ≫ 0 yields the ambipolar field with Ξ 0 = 0 and Ξ 1 = 1. In the thin transition zone |q| < ∼ λ D n i , the AAA is an approximation which is correct up to errors ∼ (λ D /l) 2 .

Once the AAA is implemented, the rest follows quickly: The ion equation of motion is also transformed into charge coordinates,

- Ψ 2 i n 2 i ∂n i ∂q = e m i E + 1 λ Ψ 2 i n 2 i . (9)
Inserting the AAA field expression (7) and defining the functions of q and n i ,

L (q, n i ) = T e m i Ξ 1 eq √ ǫ 0 T e n i - Ψ 2 i n 2 i , (10)
R (q, n i ) = 1 λ Ψ 2 i n 2 i - e m i Ξ 0 eq √ ǫ 0 T e n i T e n i ǫ 0 , (11)
one arrives at the sheath equation, an ordinary differential equation for the ion density n i as a function of the spatial coordinate q.

L (q, n i )

∂n i ∂q = R (q, n i ) . ( 12
)
The sheath equation ( 12) resembles the presheath model ( 4) but overcomes its deficiency:

Its singularity is removable (see the topology in fig. 3): Exactly one physical solution passes through the intersection (q * , n * i ) of the curves L (q, n i ) = 0 and R (q, n i ) = 0. Once the ion density n i (q) is found, the electric field E(q) follows from (7). Also, the transformation back into the x-coordinates can be carried out. Then the potential Φ(x) and the electron density n e (x) can be calculated. Fig. 4 shows the results for the parameters of the original model.

The agreement is excellent, the relative deviation is in the percentage range. (It is a small fraction of (λ D /l) 2 which is between 0.04 to 0.1 for all solutions.)

The sheath equation (12) achieves what was demanded above: It delivers a simultanous decription of both the ambipolar plasma (for q ≫ q * ) and the unipolar sheath (for q ≪ q * ), and it was derived without adopting the assumption that the Debye length is small compared to the mean free path. It is interesting to note, however, that the location of the removable singularity which links the branches depends solely on the ratio of those two parameters.

For a concise notation, the functional inverse of Ξ 0 (ξ

)/Ξ 1 (ξ) is defined Ξ(•) = (Ξ 0 /Ξ 1 ) -1 (•), its functional composition with Ξ 1 (ξ) is written V (•) = Ξ 1 (Ξ(•)) ≡ Ξ 1 ((Ξ 0 /Ξ 1 ) -1 (•)).
Denoting by an asterisk that λ * D and λ * are evaluated at the singularity itself, the reduced location ξ * = eq * / ǫ 0 T e n * i and the corresponding ion speed |v * i | = Ψ i /n * i are given as

ξ * = Ξ λ * D λ * , (13)
|v * i | = V λ * D λ * T e m i . (14)
The curves are displayed in fig. 5. For small values of λ * D /λ * , the singularity is located in the quasineutral zone (ξ * ≫ 0) and its speed |v * i | is equal to the Bohm speed v B = T e /m i . With growing λ * D /λ * , the singularity approaches the electron edge (ξ * ≈ 0) and |v * i | decreases. In the collisional regime, finally, the singularity is located deep in the depletion zone (ξ * ≪ 0) and |v * i | becomes very small. Based on these observations, and on the fact that it directly replaces the essential singularity of the presheath model ( 4), it is proposed to identify the removable inner singularity (q * , n * i ) with the sought-after 'collisionally modified Bohm point' and to call its physical characterization -in particular the ion speed condition eq. ( 14)a 'collisionally modified Bohm criterion'.

Clearly, it is of interest to compare the new criterion with the results obtained by others. Fig. 5 displays the expressions (1+λ D /λ) -1/2 (Godyak [6]), (1+(λ D /λ) 2/3 ) -1/2 (Valentini [8]), and (1 + 12.9 (λ D /λ) 0.8 ) -1/2 (Chen [9]), all proposed as the ratio of a collisionally modified Bohm speed to the standard value. (Formulas converted to the notation of this manuscript.)

All curves agree qualitatively with the curve V (λ * D /λ * ) presented here; they state that the ion speed at the edge of a collisional sheath lies below the Bohm speed, in a way that solely depends on the ratio of the Debye length λ D to the mean free path λ. There are, however, considerable quantitative differences. Undoubtedly these differences reflect a certain measure of arbitrariness in the various approaches concerning the 'exact location' of the sheath edge (and possibly also uncertainties concerning the exact definition of λ D and λ).

To what an extent, now, is the sheath criterion proposed in this manuscript not arbitrary, and how does it cope with the cited criticism of Franklin and Riemann? First, the proposed criterion is unambiguous: It relies on the characterization of a well-defined singularity within a well-defined description (10) -( 12) of the sheath. Second, the new criterion is consistent:

The description (10) - (12) was derived from the standard model ( 1) -( 3) by invoking the advanced algebraic approximation. This in turn rests solely on the assumption λ D ≪ n i / ∂n i ∂x which is safely applicable for all ratios λ D /λ. (Attested by the nearly perfect match of the density curves in figs. 1 and4.) And third, it is entirely reasonable to call (14) a 'collisionally modified Bohm criterion': As the original Bohm criterion, it marks the transition from the ambipolar diffusion dynamics of the plasma bulk to the unipolar motion in the sheath, and it reduces to its prototype in the collisionless limit. Moreover, the original Bohm point of the presheath equation ( 4) and the collisionally modified Bohm point of the sheath equation (14) proposed here are in direct mathematical analogy. In summary, this manuscript claims that there is a collisionally modified Bohm criterion:

Exactly the opposite of what Riemann and Franklin concluded. However, the difference in opinion is only about words, not about physics: Riemann's and Franklin's 'Bohm criterion' and the 'Bohm criterion' of this work rely on different assumptions and are different entities.

They are complementary to each other, not contradictory: In particular, Riemann's sheath analysis and Bohm criterion become exact in the limit λ D /λ → 0 but are not applicable for large degrees of collisionality. In contrast, the sheath model of this work relies (as the advanced algebraic approximation itself) on the condition λ D ≪ n i / ∂n i ∂x , an assumption which is never 'exact' but always 'reasonable', i.e, fulfilled to a satisfactory degree for all ratios λ D /λ.

The same holds for the modified Bohm criterion proposed in this work.

It remains to mention that the approach presented in this manuscript can be extended to RF modulated sheaths. Several investigations which have verified the presented ideas or used them for practical work are already completed [16][17][18][19][20][21][22][23][24], others -particularly in the field of microplasmas -are being conducted. Of all the possible solutions indicated by the stream, only one is physical, i.e., monotonous and defined for all q (thick curve).

  Bohm himself studied the situation on the scale x ∼ λ D . Alternatively, one can focus on the scale x ∼ λ and employ the 'plasma approximation': Substituting quasineutrality (n e = n i ) for Poisson's equation and combining it with Boltzmann's relation yields E =this 'ambipolar field' into the equation of motion and utilizing flux conservation leads to the 'preseath model' (here formulated in n i instead of the ususal v i ):

1 FIG. 2 :

 12 FIG.2:The analytical switch functions Ξ 0 (ξ) and Ξ 1 (ξ) of the advanced algebraic approximation in dependence of their argument ξ.

3 FIG. 3 :

 33 FIG.3: Phase portrait of the sheath equation at p = 10 Pa. The topology is structured by the dashed curve L(q, n i ) = 0 and the dotdashed curve R(q, n i ) = 0. Their intersection is a removable inner singularity of saddle point type which is interpreted as the collisionaly modified Bohm point.

3 FIG. 4 :

 34 FIG.4: Density of the electrons (thin) and ions (thick) from the sheath equation, for the same parameters as the original sheath model (fig.1). The deviation is less than the width of the lines.The vertical lines indicate the electron edge s. The solid points denote the generalized Bohm point defined in this manuscript, i.e., the location of inner singularity of the sheath equation.

13

 13

6 FIG. 5 :

 65 FIG. 5: The ratio V of the collisionally modified Bohm speed to the standard Bohm speed v B as a function of λ * D /λ * (solid), in comparison with the formulas presented by Godyak [6] (dashed), Valentini [8] (dotted), and Chen [9] (dotdashed). The inset shows the location Ξ of the collisionally modified Bohm point as a function of λ * D /λ * .

Acknowledgment. -The author gratefully acknowledges support by the Deutsche Forschungsgemeinschaft via FOR 1023 and SFB-TR 87.


Frédéric Roig

Eric Dantras
email: dantras@cict.fr

Jany Dandurand

Colette Lacabanne

Influence of hydrogen bonds on glass transition and dielectric relaxations of cellulose

Keywords: 77, 22, Gm, 77, 84, Jd

The molecular dynamic in hydrated cellulose has been investigated by a combination of thermal analyses and dielectric spectroscopy. Differential Scanning Calorimetry (DSC) shows the dependence upon hydration of the glass transition temperature T g . A physical ageing phenomenon has been observed. At the molecular scale, bound water is hydrogen bonded to polar sites of cellulose macromolecules. At the macroscopic scale, water molecules play the role of a plasticizer for cellulose lowering its T g . Dynamic Dielectric Spectroscopy (DDS) combined with Thermostimulated Currents (TSC) have allowed us to follow more localized molecular mobility. The β relaxation mode is characterized by activation entropies that vanish for higher water contents indicating molecular mobility localization. It is plasticized by water like the glass transition. This analogy is explained by a common origin of both mechanisms: the mobility of the cellulose backbone. The evolution of the γ mode upon hydration follows an anti-compensation law. Water acts as an anti-plasticizer in a hydrogen bonded network.

Introduction

Cellulose is probably the most studied biopolymer. Increasing CO 2 stimulate photosynthesis of chlorophyll plants and potentially allow plants to increase production biomass in the order of 20% [1].

The major structural component of the plant cell and the most important substance produced by plants is cellulose [START_REF] Wagner | Chemistry Ultrastructure Reactions[END_REF]. The large number of studies published on its properties is explained by important interests to optimize the natural resource consumption [START_REF] Gani | [END_REF]4], prevention of fire [5], bone regeneration [6], specific cellulose membranes [7,8] or decrease the cost of cellulose extraction process by wood genetic modification [9] in paper industry for example.

Thermal stability has been characterized by thermogravimetry analysis (TGA) [10] and physical structure by differential scanning calorimetry (DSC) [11]. In the literature, thermogravimetry analyses of cellulose point out two weight losses; the first one in the range of 40 to 80 °C attributed to the water evaporation [11][12][13] and the second one in the temperature range of 200 to 400°C, to polymer degradation. During this decomposition, the principal products are water and carbon dioxide [5,11,12,14,15].

The glass transition of cellulose has been extensively studied [11,[16][17][18], generally by differential scanning calorimetry [19]. The determination of the glass transition temperature is not easy due to the presence of water [16] and to the weak fraction of amorphous phase [START_REF] Hatakeyama | Structure and properties of the amorphous region of cellulose Cellulose Structural And Functional Aspects[END_REF]. Similar problem has been already noted for starch-based materials and various factors such as sample preparation, methods of adding water, effect of sample mass or heating rate influence have to be carefully considered [START_REF] Yu | [END_REF]. According to Hatakeyama [START_REF] Hatakeyama | Structure and properties of the amorphous region of cellulose Cellulose Structural And Functional Aspects[END_REF], physico-chemical properties of cellulose are mainly dependent upon amorphous regions. By x-ray diffractometry (XRD) the amorphous content of various kind of cellulose is measured. It reaches 30% for natural cellulose extracted from cotton linter [START_REF] Hatakeyama | Structure and properties of the amorphous region of cellulose Cellulose Structural And Functional Aspects[END_REF].

Batzer and Kreibich [16] determine Tg by DSC but they have small changes in heat capacity.

Szczesniak [11] has also small variation of heat capacity. Vittadini et al. [17] do not observe T g but molecular transition is found by nuclear magnetic resonance (NMR). Dielectric dynamic spectroscopy (DDS) has been already used on polysaccharides like cellulose and its derivatives [22][23][24][25][26][27][28][29]. Four dielectric relaxations modes are usually observed γ, β, β wet and α, in the order of increasing temperatures [22][23][24][25][26][27][28][29]. The dielectric properties of cellulose have been also investigated by combined DDS/TSC [28]. In the low temperature range (T < 0 °C) where the molecular mobility is localized two modes are extracted; i.e. γ and β also observed in mechanical spectroscopy [30][31]. The γ mode has been associated with side groups (-CH 2 OH and -OH) mobility [22][23][24][25][26][27][28] and the β mode with localized cooperative molecular mobility of the main chain [28,29]. The β wet relaxation mode is observable at room temperature and is not detected after intense drying. This mode is attributed by Einfeldt et al. [22][23][24][25][26][27] to the reorientation of polymer-water complexes formed at the polymer-water interface. At higher temperature (T>0°C) the molecular mobility is delocalized. α relaxation mode is assigned to cooperative delocalized mobility of nanometric sequences of the main chain [28]. It is the dielectric manifestation of glass transition.

According to Einfeldt et al., the properties of polysaccharides as cellulose can be changed drastically by adding small amounts of water. Moreover the ability of cellulose to interact with water is not reversible, after consecutive dryings the water retention value of cellulose decreases drastically. This phenomenon is known as hornification [22,23,26]. The influence of hydration on relaxation modes of polysaccharides is noticeable. Einfeldt compares dry and wet celluloses. The activation energy E a of the local main chain dynamics increases contrary to the pre-exponential factor τ 0 which is reduced. The γ process is more cooperative in wet celluloses than in dry celluloses [25]. Measurements on others polysaccharides as starch and dextran confirm this results [22]. A β wet relaxation mode can be observed in all wet polysaccharides in the room temperature range (0-60°C) [22][23][24][25][26][27]. This mode is not observed anymore after drying the sample [22][23][24][25][26][27]. The hydration influence of polysaccharides on their dielectric properties can be interpreted by several hypotheses: one the one hand, the OH-groups of glucosidic units form hydrogen bonds with water molecules. These bound water molecules increase the dipolar moment and the moment of inertia of movable groups. Moreover the stiffness of the polymer chain increases because water molecules bonded produce a bridge along the chain and between two adjacent chains. On the other hand, the macroscopic flexibility of cellulosic fibrous sample becomes more important thanks to the swelling power of water [START_REF] Wagner | Chemistry Ultrastructure Reactions[END_REF]22,25].

In spite of the abundant literature, some points as the complex interaction of water with cellulose macromolecule remains particularly obscure and multiple interpretations exist [32,33,34].

We have already shown that the fruitful synergy of DDS/TSC combination allows us to investigate the changes within the molecular dynamics of cellulose caused by water influences [28,29]. These results stimulate us to study the issue of cellulose glass transition with different water contents using the classical techniques of DSC. A complete dielectric study will be performed and particular attention will be paid to specific interactions between polymer and water.

Experimental

Materials

The cellulose provided by Sigma-Aldrich is a white microcrystalline powder (20 µm) extracted from cotton linters, with a number-average molecular weight between 36000 and 40000. The degree of polymerization is between 222 and 246. The initial water content has been determined by thermogravimetry experiments near 5%. Cotton fibres extracted from a natural plant have been also investigated (purity 95%).

Samples with different hydration levels were prepared. The 0% hydration was obtained according to the following procedure: sample was placed on P 2 O 5 during 72 h. Constant humidity level was reached by storage at 25 °C over saturated salt solutions during 72 h. Saturated salt solutions of MgCl 2 , 6H 2 O and NaCl allow us to obtain 3% and 7% of hydration respectively. In the case of NaCl solution, for a double storage time, 13% was reached.

The water content was determined from the weight loss after drying. The hydration level h(%), was defined as follows:

( )

100 % × = cellulose water m m h (1)
Where m water is the water mass and m cellulose is the mass of cellulose sample after drying.

Methods

Thermogravimetry

Thermogravimetry analyses were performed with a TGA Q50 Thermal Analyser in open pans under nitrogen flow from room temperature to 600 °C at 10°C/min. Sampling of cellulose was twofold: powder or pellets. Activation energy of cellulose degradation is analysed using Broido method [35]:

( ) ( ) ( ) RT E K ∆ - = - - ln 1 ln ln α (2) ∞ - - = m m m m 0 0 α (3)
Where m 0 is the initial mass and m ∞ is the final mass of the sample.

Standard differential scanning calorimetry

The DSC thermograms were recorded with a DSC 2920 TA Instrument. The samples in pellets gave better results than powder (associated with signal-to-noise ratio). Cellulose powder was pressed under controlled pressure of about 75 MPa to form pellets. Samples were sealed into aluminium pans, and empty pans were used as reference. Empty pans were also sealed. Experiments were performed in the -20 to 180°C temperature range with a heating rate 20°C/min (and a cooling rate q c = 20°C/min).

This heating rate has been chosen to perform the thermal signal [36].

Dynamic dielectric spectroscopy (DDS)

A BDS 4000 Novocontrol broadband dielectric spectrometer system was used to obtain the dielectric relaxation map in a wide temperature and frequency range. The sample was inserted between two plan parallel electrodes. It was stimulated by a sinusoidal electrical field at a given temperature and the complex impedance is recorded during frequency scan. The measurements were carried out in the frequency range of 10 -1 -10 6 Hz from -150°C to 200°C by steps of 5°C. The complex dielectric permittivity ε * (ω) was deduced thanks to complex impedance by the relation:

( ) " ' ) ( * 1 * 0 ε ε ω ω ω ε i Z C i - = ⋅ = (4)
With C 0 is the capacitance in vacuum, ω the pulsation (ω=2πf), Z * the complex impedance. Relaxation modes are described by the Havriliak-Negami function [37]:

( ) ( ) [ ] HN HN HN i β α ωτ ε ε ω ε + ∆ + = ∞ 1 * (5)
With τ HN is the relaxation time of the Havriliak-Negami (HN) equation, α HN is the width and β HN the asymmetry of the mode.

In cellulose, the temperature dependence of relaxation times follows an Arrhenius-Eyring equation:

∆ ⋅ ∆ - ⋅ ⋅ = RT H R S T k h T B exp exp ) ( τ (6)
Where k B is the Boltzmann's constant, R is the gas constant; h is the Planck's constant, ∆H the activation enthalpy and ∆S the activation entropy.

Thermo stimulated currents (TSC)

TSC measurements were carried out with a home-made equipment previously described [19]. Pellets (φ = 0.8 cm) were inserted between two plate steel electrodes. Sample cell was filled with dry helium.

For recording complex thermograms, the sample was polarized by an electrostatic field E p = 1.2 kV/mm during t p = 2 min over a temperature range from the polarization temperature T p down to the freezing temperature T 0 . Then the field was turned off and the depolarization current was recorded with a constant heating rate (q h = + 7 °C.min -1 ), the equivalent frequency of the TSC spectrum was f eq 10 -2 -10 -3 Hz. Elementary TSC thermograms were obtained with a poling window of 5°C. Then the field was removed and the sample cooled down to a temperature T cc = T p -40°C. The depolarization current was normalized to be homogeneous with dipolar conductivity σ and it was recorded with a constant heating rate q h .The series of elementary thermograms was recorded by shifting the poling window by 5°C towards higher temperature.

Results and discussion

Thermal stability

TGA experiments have been carried out in order to define the stability and the influence of sampling (powder or pellet). Figure 1 shows two thermograms of powder and pellet respectively. For powder two weight losses are observed, i.e. 60°C and 343°C. First one is associated with water according to infra-red spectrometer and second one with the cellulose degradation. These two events are consistent with literature [4,5,[11][12][13][38][39][40]. For pellets, a slight shift of these two weight losses toward higher temperatures is observed, i.e. 82°C and 355°C. This shift is probably due to the increase of hydrogen bonds interactions associated with the decrease of length between hydrophilic groups.

Broido method [35] has been used to determine activation energy of degradation. Results are shown in Figure 2 and reported in Table 1. Three sampling of cellulose (i.e. powder, pellet and fibres) and three pressures are analysed. The activation energy of pellets, near 238 kJ.mol -1 , is higher than powder and fibres. This value is consistent with Varhegyi et al data [41]. In the case of pellets, we do not observe influence of pressure on this physical parameter. Analogous behaviour has been found in tablets of microcrystalline cellulose [42]. Consequently, pellets can be used for dielectric analysis. The stability temperature range as defined by TGA reaches 200°C.

Thermal transitions

DSC thermograms of cellulose are reported in Figure 3. Four consecutive heating cycles have been carried out to obtain a repeatable result in dehydrated cellulose. A large endothermic peak associated with water loss during the heating [11,43] is observed for the first thermogram. For the following thermograms the heat capacity step characteristic of the glass transition is pointed out. The T g of dehydrated cellulose is found at 84 ± 1°C with a ∆C p = 0.30 ± 0.01 J/(g.°C). The ∆C p value is comparable with values measured for synthetics polymers [START_REF] Wunderlich | Thermal Analysis Of Polymeric Materials[END_REF].

Figure 4 reports an endothermic peak superimposed to T g appearing when the sample was heated at 60°C for different duration from 15 min to 3 h. This weak phenomenon is reversible so that it looks like physical ageing [START_REF] Struik | Physical Aging in Amorphous Polymers and Other Materials[END_REF][START_REF] Samouillan | [END_REF]. Its kinetics must be fast since it is independent from annealing time contrarily with usual physical ageing. It is interesting to report here that physical ageing has been observed by thermo stimulated currents in amorphous ethyl cellulose [47]. At this stage, the phenomenon observed in cellulose may be attributed to the disruption of hydrogen bonds in the amorphous phase.

Figure 5 shows heat flow measurements of cellulose with water content from 0% to 5%. The glass transition is strongly dependent upon hydration. T g and ∆C p of cellulose decreases when the hydration level increases. A strong plasticization effect is observed and confirms previous results [11,16,18].

The properties of biopolymers are dependent on their hydration due to hydrophilic groups [18]. Authors differentiate three kind of water in biopolymers: free water, freezing-bound water and non-freezing bound water [11,42,48,49].

This classical thermal analysis reveals that water interacts with the polymer; i.e. free water is dismissed. The following dielectric analysis allows us to shed some light on this behaviour.

Isothermal dynamic dielectric response

The complex dielectric permittivity is explored by DDS. In order to investigate the water influence on the dielectric response, experiments are carried out with different hydration levels.

Figure 6 shows the relaxation map of 7% hydration cellulose. These three-dimensional diagrams point out four dielectric events: at low temperature two broad secondary relaxation modes called γ and β, at higher temperature a wide quasi-isothermal relaxation and conductivity phenomenon labelled β' and σ respectively. The α relaxation mode associated with the dielectric manifestation of glass transition is hidden by this conduction phenomenon. The σ electrical conduction manifestation has been already associated with charge transport [26] or Maxwell-Wagner-Sillars effect [28].

The tan (δ) loss factor of cellulose with different water content; i.e. 0% to 7% is shown in Figure 7. In the low temperature region, γ relaxation mode is pointed out for each hydration level contrary to β relaxation mode which is only visible for 3 and 7%. In consequence for 5% water contents and near -50°C, the tan (δ) loss factor curve is below 3%. At higher temperatures (>80°C) a similar behaviour is observed due to the influence of α relaxation mode hidden by conduction phenomenon. From -50°C to 80°C the additional dielectric contribution labelled β', appears for hydrated samples. This peak is also observed on the three-dimensional map of ε'' (Figure 6). The strength of this mode increases with the hydration level. Such a β' relaxation mode, has been already reported in the literature [42,50,51] and it has been associated with conformational changes of confined water. Banys et al. [50] confirms this molecular origin in porous materials.

The relaxation times τ have been extracted from dielectric data and reported in the Arrhenius diagram of Figure 8. The relaxation times of the two secondary relaxation modes are well fitted by an Arrhenius law but their behaviour is dependent upon hydration level. For the γ mode, the series of relaxation times corresponding to various hydration levels converges at low temperature and follows an "anti-compensation" law:

- ∆ ⋅ = C C T T R H T 1 1 exp ) ( τ τ (7)
Where T c = -130°C and τ c = 2.4 s are the temperature and the time of the convergence point respectively. This τ c value is quite high for a low temperatures relaxation. As comparison, it is near τ c value of the lignin main relaxation [52] It is interesting to note the value of ∆H= 40 kJ/mol is comparable with the one of hydrogen bond. The activation parameters are reported in Figures 9 and10. These results are coherent with published data from mechanical and dielectric studies [28][29][30][31]49]. The activation energy of the γ mode increases with the hydration level whereas it decreases for the β mode. The increase of the activation energy indicates the decrease of its molecular mobility [23,28,53]. The hydroxyl and hydroxymethyl groups form water-polymer complexes. Hydrogen bonds decrease the water-polymer complex mobility. Hydration plays the role of an anti-plasticizer for the γ mode. The decrease of the activation energy with hydration for the β mode indicates that its molecular mobility increases. The water molecules interfere with existing hydrogen bonds, increasing the free volume between the polymer chains and reducing the hindrance for the local chain mobility. Hydration plays the role of a plasticizer for the β mode. It can involve in the same behaviour for the β mode of epoxy assemblies associated with hydroxyether entities. Chevalier et al. [54] observe also the low temperature convergence point and the decrease of the activation energy with the increasing of hydration level due to hydrothermal ageing. The evolution of mobility is linked with the physical environment modification. These dynamic dielectric data have been completed towards lower frequencies by thermostimulated current experiments.

Correlation with non isothermal dielectric response

The normalized depolarization current σ is plotted versus temperature in Figure 11, for different water content. In the low temperature range and for 0% hydration two relaxation modes are observed. The γ relaxation mode assigned to the side groups' reorientation is located at -135 °C. The β relaxation mode attributed to localized movements of short sequences of the main chain of cellulose is pointed out at -65 °C. Contrary to the β relaxation mode, the maximum of the γ relaxation mode increases with the hydration level from -135 °C at 0% to -120 °C at 7%. We observe a merging of these two relaxation modes upon increasing hydration level. These results are coherent with SDD; i.e. hydration plays the role of an anti-plasticizer for the γ mode and a plasticizer for the β mode. The antiplasticization phenomenon due to hydration is also observed in micromechanical spectroscopy as shown by Lamure et al. from studies of the main retardation mode associated with glass transition of proteoglycan subunits [55]. Fractional polarizations allow us to extract activation parameters for different hydration state. A series of elementary thermograms is obtained and reported in Figure 12 for 0% hydration.

Results of the two dielectric methods are reported for each hydration state. The β DDS relaxation mode is detectable for 3% and 7% hydration level only. According to Montès et al. [30,31,56], the β relaxation mode has been associated with cellulose moisture content; it is not always detectable by DDS. In Fig. 13, Arrhenius diagram of DDS and TSC data is shown for 3% of water.

The relaxation times obtained by DDS are extrapolated until TSC data; a good correlation is performed. The correlation between γ and β modes in DDS and TSC allows us to confirm the molecular origin of these two modes. Then, the β mode can be follow for each state of hydration by TSC.

The activation enthalpy versus temperature for γ and β relaxation modes are reported in The influence of hydration level on the activation enthalpies versus temperature for the elementary process of γ and β relaxation modes are reported in Figure 15. The enthalpy range and ∆H max values decrease when the hydration level increases. The cooperative behaviour of the β relaxation is strongly dependent upon hydration. For water contents between 0% and 2.5%, the β relaxation is characterized by a cooperative behaviour; otherwise for hydration level between 2.5% and 7%, a non cooperative behaviour is exhibited. This evolution has been assigned to an increase of the hydrogen bond density leading to a more localized cooperative mobility.

Conclusion

The combination of both thermal analyses and dielectric spectroscopy has been applied to a thorough investigation of the molecular dynamic in hydrated cellulose.

The thermal stability of the polymer and the initial water content is assessed using TGA.

Thermal behaviour of cellulose is dependent of morphology of samples. The activation energy of degradation for pellets is higher than for the others forms (fibres, powder). Contrarily, pressure does not affect the degradation behaviour of pellets.

Original data have been extracted from DSC showing the influence of water on the viscoelastic behaviour of cellulose through the analysis of the dependence upon hydration of the glass transition temperature T g . At the molecular scale, bound water is hydrogen bonded to polar sites of cellulose macromolecules. At the macroscopic scale, water molecules play the role of a plasticizer for cellulose lowering its T g . DDS combined with TSC has allowed us to follow more localized molecular mobility occurring in the glassy state as a function of hydration. According to the literature, the origin of the β relaxation mode is localized movements of the β 1-4 glycosidic bonds of cellulose. This mode is characterized by activation entropies that vanish for higher water contents indicating a localization of the molecular mobility. It is plasticized by water like the glass transition. This analogy is explained by a common origin of both mechanisms: the mobility of the cellulose backbone. In fact, they differentiate only by the scale of the mobile units.

From previous works, the γ relaxation mode has been attributed to the molecular mobility of the side groups of glycosidic rings. The behaviour of this γ mode upon hydration is opposite to the previous ones. The series of relaxation times corresponding to various hydration levels follow an anticompensation law with an activation enthalpy comparable with the one of hydrogen bond. Water acts as an anti-plasticizer in a hydrogen bonded network involving -CH 2 OH and -OH groups with water molecules. It is interesting to note here that such an arrangement is common to biopolymers containing glycosidic rings.

Figure 14

 14 Figure 14 with the Starkweather line [57] associated with the null activation entropy. The experimental points related to the γ mode, lie close to the Starkweather line indicating that the localized mobility of the side groups is not cooperative. Contrarily, activation enthalpies of the β mode

Figure 1 :

 1 Figure 1: Thermogravimetry thermograms for cellulose powder and pellet, from room temperature to 550°C at 10°C/min.

Figure 2 :

 2 Figure 2: Broido's representation for different samples of cellulose.

Figure 3 :

 3 Figure 3: Heat flow thermograms of cellulose pellets for four successive scans from 0°C to 180°C at 20°C/min.

Figure 4 :

 4 Figure 4: DSC thermograms of cellulose after different annealing times from 20°C to 180°C at 20°C/min.

Figure 5 :

 5 Figure 5: DSC thermograms of hydrated cellulose from 30°C to 120°C, recorded at 20°C/min. The hydration level h is indicated for each thermogram.

Figure 6 :

 6 Figure 6: Imaginary part of dielectric permittivity versus temperature (-160 °C to 200°C) and frequency (10 -1 Hz to 10 6 Hz), for cellulose pellets with 7% water content.

Figure 7 :Figure 8 :

 78 Figure 7: Tan (δ) dielectric loss factor versus temperature (-160 to 200°C) at 1Hz for cellulose at different hydration levels

Figure 9 :

 9 Figure 9: Activation energy of γ and β mechanical and dielectric relaxation modes for different moisture contents for cellulose.

Figure 10 :

 10 Figure 10: Activation entropy of γ and β mechanical and dielectric relaxation modes for different moisture contents for cellulose.

Figure 11 :

 11 Figure 11: Complex TSC themogram of cellulose at four states of hydration from -160°C to -10°C at 7°C/min.

Figure 12 :

 12 Figure 12: Complex TSC themogram with their associated elementary thermograms for γ and β relaxation modes (-160°C to 0°C at 7°C/min) for cellulose in the dehydrated state

Figure 13 :Figure 14 :

 1314 Figure 13: Arrhenius diagram from combined DDS and TSC in the low temperature range (T< 0°C) for 3% hydrated cellulose

Figure 15 :

 15 Figure 15: Activation enthalpy versus temperature for the elementary thermograms constituting the low temperature relaxation modes for cellulose at four states of hydration.

Figure 1 (Figure 2 Figure 16 (

 1216 Figure 1 (Figure1.tif)

Table 1 .

 1 Activation energy of degradation for different cellulose samples Cellulose

	Fibers	120 ± 3
	Pellet (P=150 MPa)	238 ± 10
	Pellet (P=75 MPa)	237 ± 13
	Pellet (P=37.5 MPa)	232 ± 4
	Powder	225 ± 4
		Ea (kJ/mol)


M A Sbai
email: a.sbai@brgm.fr

M Azaroual
email: a.azaroual@brgm.fr

A Numerical Model for Miscible Displacement of Multi-Component Reactive Species

The study of instable miscible displacement of chemical species is an important process for many environmental and reservoir engineering applications. A new coupled saturated porous media ‡ow and reactive transport numerical model: RTAFF, has been developed for the assessment and analysis of complex feedback loops between physical transport, geochemistry, and ‡uid ‡ows in such environments.

The model accounts for changes of ‡owing phase mobility upon viscosity of individual species in the aqueous solution, and …eld permeability attributed mainly to mineral precipitation-dissolution reactions. An operator splitting framework sequentially couples a mixed-hybrid …nite element ‡ow model on quadrilateral/hexahedral grids, with a mass conservative advective-dispersive-reactive (ADR) model able to handle a full range of homogeneous and heterogeneous reactions, using a mixed kinetic equilibrium approach. The multi-component species ADR simulator work ‡ow involves a time-split procedure adapted for each single process solver, thus advection and dispersion operators are also solved independently using other specialized computer models. Example applications are presented in one, and two space dimensions, showing the e¤ectiveness of the numerical model for practical environmental applications. First results demonstrate the accuracy, e¢ ciency and robustness of our approach. Therefore, the model may prove useful for many practical applications.

INTRODUCTION

Demands to undertake modeling analysis of coupled groundwater ‡ow, solute transport, and reactive water-rock interactions, is increasing among scientists and engineers in di¤erent disciplines. The development of numerical models for coupled subsurface phenomena is growing considerably, as current literature attests [START_REF] Lichtner | Continuum formulation of multicomponent-multiphase reactive transport[END_REF][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][START_REF] Sbai | RTAFF -A new computer model for reactive transport and ‡ow in saturated media[END_REF][START_REF] Steefel | A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive ‡ow in single phase hydrothermal systems[END_REF][START_REF] Steefel | Approaches to modeling reactive transport in porous media[END_REF][START_REF] Yeh | A model for simulating transport of reactive multispecies components: model development and demonstration[END_REF]. Several numerical di¢ culties may arise when many interacting phases, like dissolved aqueous species, immobile adsorbed species, and precipitation / dissolution of minerals and gases, are taken simultaneously, and are a¤ecting each other. Concentrations of chemical species, even for those present in small amounts, must be kept positive. Mass conservation has to be veri-…ed locally for each grid cell. Non-oscillatory behavior for highly advective ‡ows should be suppressed, to accurately capture physical sharp moving fronts of, aqueous components, and mineral phases. Therefore, stringent requirements are needed on the performances of numerical techniques.

In this paper we present RTAFF, which is a new reactive transport model developed recently. In the process of designing the simulator, several issues related to (i) robustness, (ii) accuracy, and (iii) e¢ ciency are considered. Another focus is to illustrate the challenges and di¢ culties in applying the presented numerical algorithms to practical problems. Realistic applications involve quite often multi-dimensional domains, highly heterogeneous permeability …elds, and reactions spanning di¤erent time scales.

MATHEMATICAL MODEL

The mathematical model is a combination of three modules for, groundwater ‡ow, solute transport, and batch geochemistry. An overview of the governing equations, with the underlying approximations taken are brie ‡y presented in the subsections below. A complete description of the model should be found elsewhere [START_REF] Sbai | RTAFF -A new computer model for reactive transport and ‡ow in saturated media[END_REF].

Groundwater ‡ow

The governing equations for three-dimensional, saturated, incompressible, groundwater ‡ow are derived from the (i) mass conservation, and (ii) generalized Darcy's laws, given as (see [START_REF] Bear | Dynamics of ‡uids in porous media[END_REF][START_REF] Sbai | Modeling three-dimensional groundwater ‡ow and transport by hexahedral …nite elements[END_REF] for more details)

S 0 @ @t + r:u = q (1) u = k [r grz] (2)
where S 0 is the speci…c storage coe¢ cient [L 1 ], is the pressure potential [L], u is the Darcy velocity vector [L/T], q is the external source and/or sink term [T 1 ], k is the permeability tensor [L 2 ], is the dynamic ‡uid viscosity [M/LT], g is the acceleration due to gravity [L 2 /T], is the ‡uid density [M/L 3 ], r = @ @x ; @ @y ; @ @z T is the del operator

[L 1 ]
, where x = (x; y; z) T is the location vector [L], the z-axis is taken vertical and positive downwards, and t is time [T].

Reactive geochemical transport

Mass conservation laws in saturated porous media, for a multi-component reactive system, yield the following nonlinear partial di¤erential equations (see [START_REF] Lichtner | Continuum formulation of multicomponent-multiphase reactive transport[END_REF][START_REF] Steefel | A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive ‡ow in single phase hydrothermal systems[END_REF][START_REF] Yeh | A model for simulating transport of reactive multispecies components: model development and demonstration[END_REF] for more details)

@( U j ) @t + r: [J Ad + J Disp ] = R K j (C 1 ; :::; C N S ) + q j j = 1; :::; N C (3)
where, N S is the total number of species, divided into N C primary basis species and N R secondary species set, e.g. the number of independent chemical reactions [START_REF] Lichtner | Continuum formulation of multicomponent-multiphase reactive transport[END_REF], is the porosity, C i is the concentration of species i [moles/L 3 ], U j is the total aqueous concentration of primary species j in solution [moles/L 3 ], J Ad = uU j and J Disp = DrU j , are the advective, and dispersive ‡uxes respectively [moles/L 2 T], R K j = dC j dt is the total reaction rate of all kinetic reactions involving species j [moles/L 3 T], D is the hydrodynamic dispersion tensor [START_REF] Bear | Dynamics of ‡uids in porous media[END_REF], including the e¤ects of molecular di¤usion coe¢ cient, longitudinal, and transverse dispersions [L 2 /T]. q j is the source and/or sink term [moles/L 3 T] representing injection or production wells, and expressed as

q j = jqjU ext j jqjU j (4)
where, jqj is the absolute value of the positive/negative injection/extraction ‡ow rate, and U ext j is an external total concentration. Chemical reactions are either su¢ ciently fast, allowing for a thermodynamic equilibrium description, or rate-limited. Examples of the N Q R equilibrium reactions, are the homogeneous aqueous phase, complexation reactions, or surface complexation reactions. Heterogeneous aqueous-mineral reactions, are part of the N K R kinetic, or rate-limited reactions. Notice that,

N R = N Q R + N K R .
Total concentrations are de…ned as:

U j = C j + N Q R X i=1 j;i+N C C i+N C j = 1; :::; N C (5)
where, j;i is the stoichiometric coe¢ cient for species i in component j. Concentrations of secondary aqueous species are determined from mass action laws:

8 < : log a i+N C = N C P j=1 j;i log a j + log K i i = 1; :::; N Q R log a i = log i + log C i i = 1; :::; N C + N Q R (6)
where, a i and i are the molal activity and activity coe¢ cient of species i respectively. We use Davies equation [START_REF] Morel | Principles and applications of aquatic chemistry[END_REF] for low ionic strength solutions, I, and given by

log i = Az 2 i p I 1 + p I 0:3I ! (7)
where, z i is the ionic charge of aqueous species i, A is a temperature dependent co-e¢ cient, and the ionic strength is I = 1 2

Naq P i=1 z 2 i C i .
Davies equation links the activity coe¢ cients with the concentrations of aqueous species. Concentrations of equilibrium mineral phases are given by: 8 > > > < > > > :

N C P j=1 j;i log a j + log K i = 0 if a i+N C > 0 N C P j=1 j;i log a j + log K i 6 0 if a i+N C = 0 (8)
It is interesting to highlight the non di¤erentiable character of the nonlinear algebraic equilibrium equations 6 and 8. When kinetics are involved, the problem becomes of coupled di¤erential algebraic nature [START_REF] Steefel | A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive ‡ow in single phase hydrothermal systems[END_REF][START_REF] Steefel | Approaches to modeling reactive transport in porous media[END_REF][START_REF] Yeh | A model for simulating transport of reactive multispecies components: model development and demonstration[END_REF].

Coupled ‡ow and reactive transport

Groundwater ‡ow is coupled with reactive transport to allow for simulation of miscible displacement processes in contaminant transport, and reservoir engineering applications. In this case, we allow a variable mixture viscosity in the form:

= F " 0 ; i j i=1;:::;N S ; i 0 i=1;:::;N S # (9)
In general the function F depends on the resident ‡uid viscosity 0 , the mass fraction of individual species in solution, i = C i M i = , where M i is the molecular weight of species i, and i = 0 is the viscosity ratio or the inverse of the mobility ratio.

Such expressions are given as a …fth order polynomial at given temperature and pressure values, or as a power quarter mixing rule [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF]. The later is traditionally used in petroleum reservoir engineering, and expressed as

0 = " N S X i=0 i 0 1=4 i # 4 (10)

THE NUMERICAL MODEL

Solving the coupled set of equations 1 and 3, subject to equations 2, 5, 6 and 8 in one step is notoriously a complex task for multi-dimensional problems. Although its feasibility has been demonstrated by several researchers (see [START_REF] Steefel | A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive ‡ow in single phase hydrothermal systems[END_REF] and the references therein) it remains of very high CPU cost and not suitable for integration of independent evolving computer codes. Therefore, in this study operator-splitting methods are preferred, because they enable (a) tight integration of di¤erent specialized computer models in the same framework, (b) breaking down the problem nonlinearity by isolating the geochemical batch computations at each time step, (c) and making it possible to use high Péclet number algorithms, like high resolution …nite volume methods [START_REF]Godunov methods -theory and applications[END_REF] (HRFV). As already stated by [START_REF] Steefel | Approaches to modeling reactive transport in porous media[END_REF], the use of shock-capturing HRFV methods may be considered as the chief advantage of operator-splitting techniques.

Discretization techniques and numerical methods

The numerical model hierarchy includes four integrated modules:

Mixed Hybrid Finite Element (MHFE) Flow Module

The ‡ow model is approximated by either quadrilateral or hexahedral elements of general shape. A mixed hybrid …nite element approximation solves for pressure heads simultaneously at cells centers and edges of the mesh. In a post-processing step the numerical ‡uxes at the mesh faces are recovered. Constant RT0 (Raviart-Thomas) scalar basis functions are used for approximating the pressure head at the cell center, while vector basis functions are associated with the edge ‡uxes. The time discretization scheme is fully implicit, and enhanced by the use of an e¢ cient sparse incomplete factorization preconditioned conjugate gradient solver [START_REF] Sbai | Modeling three-dimensional groundwater ‡ow and transport by hexahedral …nite elements[END_REF]. The model supports distributed cell-centered parameters for all material properties, full tensor quantities, and common boundary conditions for hydrogeological applications. Recently, it has been integrated in a modern pre-and post-processing environment, where several automatic mesh generator packages with mesh re…nement and general geometry …tted possibilities are also available.

High Resolution Finite Volume (HRFV) Advection Module

Adective transport is approximated by an explicit cell-centered higher order Godunov method. This is basically a two-step method. First, an approximate Riemann solver [START_REF]Godunov methods -theory and applications[END_REF] computes new conservative ‡uxes at cell-interfaces, from cell-centered concentrations at the previous time step. Next, a nonlinear reconstruction method with appropriate slope limiters is used to prevent oscillatory behavior near the sharp fronts. The threedimensional slope limiting function prevents local minima and maxima to rise below and above their respective values at the previous time level. Therefore, undershoots or/and overshoots do not arise during the solution procedure. This is of great importance for the next succeeding reaction step in the time splitting algorithm.

The numerical scheme allows for high Péclet numbers to be taken, and thus even coarse meshes exhibit only a small amount of numerical di¤usion. This however, comes at price of a Courant-Friedrich-Levy (CFL) constraint limiting time step size, which is enforced directly in the computer model.

Batch geochemical reactor Module

A new geochemistry module, developed from scratch [START_REF] Sbai | Un nouveau module géochimique de calcul d'équilibre et cinétiques mixtes[END_REF] and based on mixed equilibriumkinetic formulation, is beginning to prove e¢ cient and robust in coupled mode. The thermodynamic equilibrium batch engine uses robust constrained optimization techniques [START_REF] Meza | OPT++: An object-oriented class library for nonlinear optimization, Sandia laboratories report[END_REF], such as the interior points method [START_REF] Karpov | Modeling chemical mass transfer in geochemical processes: thermodynamic relations, conditions of equilibria, and numerical algorithms[END_REF] to minimize the Gibbs free energy function, subject to mass-balance and positive activities constraints. The model also enables simulation of biodegradation kinetics using various Monod rate law forms.

Dispersion-di¤ usion Module

Dispersion is approximated again by the MHFE module. Here we just reuse the ‡ow solver, with modi…ed tensor and scalar variables. All transport parameters, such as dispersivities, and molecular di¤usion may be given as cell-centered distributed data sets.

Time splitting algorithm

The ADR equations are solved using an adaptive time splitting algorithm, which occupies a central place in the computer code framework. This algorithm accounts for the possible failure of any module. When this happens we do back stepping with a smaller time step size. Let Adv and R be the number of sub-steps for solving the advection and reaction terms for a given reactive time step t respectively, and k the current time step level. The algorithm proceeds as follows:

1. Advection operator: Determine the optimum time step t CF L for the mesh, guarantying the CFL time step constraint, and hence adv . For 1; :::; adv solve explicitly the system:

U k;adv j U k j t= adv = r:U k j u k + q j f orj = 1; :::; N C (11)
2. Reaction operator:

For 1; :::; R solve the algebraic-di¤erential system:

U k;R j U k;adv j t= R = R K j C k;adv 1 ; :::C k;adv N C ; C k N C +1 :::; C k N S for j = 1; :::; N C ( 12
)
if convergence is not reached for a given reactor, then t = t=2, goto 1

Dispersion-di¤usion operator:

solve implicitly the system:

U k;disp j U k;R j t = r:D k rU k;disp j for j = 1; :::; N C (13)

Speciation and Equilibrium step:

The aqueous composition is equilibrated using the chemical equilibrium operator EQ exclusively:

U k+1 j = EQ C k;disp
We point out that looping over the grid cells is internal to the advection routine, and external to the reaction step. This is enforced to enhance data locality, such that little code restructuring will be needed in distributed computing tasks. Porting the presented algorithms to parallel computer platforms is a future extension.

DEMONSTRATION EXAMPLES

Since model development begun, in late 2002, several test problems were designed for veri…cations purposes. Some simple one-dimensional models, as the …rst example shown below, were compared with the results of other available models, to establish the accuracy, and eventually the e¢ ciency of the implemented algorithms [START_REF] Sbai | Un nouveau module géochimique de calcul d'équilibre et cinétiques mixtes[END_REF]. Based on this …rst knowledge, it is trusted that the computer model, in its present state, is able to handle complex coupled mechanistic processes, of ‡uid ‡ow, and reactive geochemistry, in multi-dimensional heterogeneous domains. Many ongoing research e¤orts, under several integrated projects, will be using the developed model as a tool for investigating environmental, and groundwater pollution problems.

One-dimensional veri…cation example

This benchmark problem was designed to assess the model capabilities to account simultaneously for combined aqueous complexation, and dissolution / precipitation reactions, including oxidation / reduction reactions. First, a comparative study with PHREEQC computer model [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] is performed. A numerical grid convergence study follows.

A rich sulfuric acid solution corresponding to an initial pH of 4 is continuously injected at the upstream boundary of a soil column. The problem domain is a 5 dm length, homogeneous, and composed mainly from gypsum, and calcite minerals. Initial and boundary conditions, as well as other simulation parameters are given in table 1. 1 Initial and boundary conditions for example 1.

Component Initial conditions (mol/l) Boundary conditions (mol/l) S(-II) 0 1E-3 S(+VI) 0 0 C(+IV) 1E-4 5E-6 Ca 1E-4 5E-6 pH 13 4 pE 4 4 Calcite 0 - Anydrite 0 - Gypse 0:1 - Table
A pore water velocity of 1:6E-6 m/s is maintained constant along the simulation run, such that the ‡ow …eld is not computed. This problem involves N C = 6 number of chemical components, and N R = N Q R = 13 independent chemical reactions. In both models, thermodynamic data is taken unmodi…ed from the LLNL database …le distributed with 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].

Numerical simulations were run on uniform grids of 500 and 250 cells, for the two models PHREEQC and RTAFF respectively. Other simulation parameters of interest are: L = 0:0067m, = 0:32, t f inal = 84000 seconds, and CFL = 0:5. Simulations showed a sequence of precipitation and dissolution of calcite. As shown in Figure 1B this zone is very narrow, and is represented by the pic concentration. This is a classical result of the strong pH increase downward the ‡ow gradient. The spatial pH zone favorable to a wide precipitation area of calcite is restricted. There is a very good agreement between the results obtained with the two models. An interesting fact to mention however, is that the 250 cells simulation was intractable with PHREEQC model.

Within RTAFF code, this simulation was run in 421 reactive time steps (the …rst corresponding to aqueous speciation), since the operator-splitting algorithm has computed a t CF L = 200 seconds for all time steps. A total CPU time of 14min20s was needed on a mobile AMD Athlon XP computer, of 1:4GHz, and with 512MB RAM. and pH, at 84000 seconds.

Acid gas migration

Acid gas injection in deep saline aquifers, and hydrocarbon reservoirs is presented as a promising way to reduce greenhouse gas emissions into the atmosphere. Geological sequestration options of carbon dioxide, and other greenhouse gases, are therefore studied by numerical models. This demonstration examples aims to show that the developed model, could be useful in understanding the involved ‡uid-gas-rock interactions.

A two-dimensional [0,100] x [0,100] meters domain is divided into a 100x100 uniform grid. A stochastic permeability distribution is obtained with the HYDRO_GEN [START_REF] Bellin | HYDRO_GEN: A new random …eld generator for correlated properties[END_REF] …eld generator, with an average value of 2m/day, and a correlation length of 5m, that is the size of 5 cells. The initial geochemical composition, as well as for the recharging groundwater from the left boundary, are identical to the ones described in table 1. Other relevant parameters for this problem are: D 0 = 10 9 m 2 /s, = 0:4, t f inal = 50 days, and CFL = 0:5. The e¤ects of macrodispersion are deliberately eliminated, creating highly advective ‡ow in permeable zones, and di¤usion controlled transport in clogged pores.

Figures 2 and3 illustrate the role of …eld heterogeneities in triggering displacement instabilities. Fingered displacement of reactive species occurs with di¤erent and contrasted mobilities. Even when …eld dispersion is neglected, the concentrations of some species are not anymore front controlled, as the intuition would suggest.

By contract to the one dimensional example, small gray spots in …gure 2B indicate disparate zones of spatial calcite distribution. These zones are non-uniformly distributed. Other simulations scenarios conducted by including dispersion, indicate that they grow up with increasing dispersivities and mixing. This …rst result would suggest that sequestration in mineral form predicted by lower dimensional models would be not su¢ cient enough for other spatial scales. Another remarkable result is that, another series of simulations have showed the dominant role of …eld permeability in controlling the dissolved species displacement. By increasing CO 2 viscosity ratio to a factor of 40, which is representative for CO 2 disposal conditions at depths below 800 meters (see [START_REF] Garcia | Flow instabilities during injection of CO 2 into saline aquifers[END_REF] and references therein), we notice negligible changes in the distribution of all species.

A typical simulation was run in 1417 reactive time steps, 217 time step cuts were necessary to capture the fastest fronts. A minimum time step of 9:6E-4 days was reached, and a total CPU time of 1h35min17s was needed, among which 1h0min35s CPU was devoted to the reactive time steps (on the same platform as in the …rst example).

The results showed here are so far preliminary. Indeed, the possible improvements may include, a kinetic parameterization of the slow mineral reactions, taking into account a more detailed speciation, and perform sensitivity analysis with respect to other key transport parameters.

CONCLUSIONS

In this study, a new numerical model has been developed for an accurate description of mechanistic coupled processes of, ‡uid ‡ows, solute transport, and chemical reactions. The computer code is made extensible through e¢ cient and ‡exible operator time splitting techniques, integrating several specialized modules, and adapted internally for each single operator. These models use di¤erent approximation techniques including, the HRFV based advection, and the MHFE dispersion-di¤usion. Numerical experiments, and other code comparison studies, have shown the e¢ ciency of the proposed approaches. Sharp aqueous and mineral fronts are accurately captured, and species concentrations are stabilized to remain positive during the solution procedure. Theses results suggests that the numerical model, may be applied to several groundwater pollution, and environmental assessment problems.

1 ;

 1 :::C k;disp N C ; C k+1 N C +1 :::; C k+1 N S for j = 1; :::; N C

Figure 1 .

 1 Figure 1. Comparison between PHREEQC results (solid lines) and RTAFF results (dashed lines), for (A) some carbonate aqueous species, (B) calcite amount (in moles)and pH, at 84000 seconds.

Figure 2 .H2SFigure 3 .

 23 Figure 2. pH and calcite mineral fronts after 20 days. (A) (B)

ACKNOWLEDGMENTS

This work is supported by the BRGM research directorate in the framework of groundwater modeling and code couplings project.


THE ECORS SYSTEM: A MOBILITY DECISION-MAKING TOOL BASED ON EARTH OBSERVATION DATA

BACKGROUND AND OBJECTIVES

In operational contexts where forces have to be deployed on unknown areas, the capacity of terrestrial vehicles to move both on and off road can be determinant. Modeling such capacity remotely, without any in situ observations, constitutes the objective of the ECORS project. Scientists and professionals of the geotechnical, remote sensing and decision-making system domains have joined their forces to develop a prototype which could help military forces in decisions related to mobility. This prototype intends to characterize operations theatres from Earth observation data, and to produce geophysical and environmental maps as well as decision aid documents. The system architecture allows a fast and an automatic analysis that can be easily implemented in cartographic representations.

APPROACH & METHODS

The parameters affecting the mobility of an ATV (All Terrain Vehicle) in soft ground interact in complex ways and depend on vehicle performance, the pilot's actions, the meteorology, the geologic parameters.... The purpose of ECORS is to develop a decision-making tool based on such parameters to estimate the vehicle performance. For different kinds of vehicles, decision maps are oriented from a strategic, operative and tactical point of view. The main data sources authorized are those coming from remote sensors, including remotesensing data, existing databases, i.e., Topobase Defense, HWSD, NCAR DB, etc or opensource data. The elementary layers produced from these data and used in the calculations are: the tracks and roads, the terrain slope, the landuse, the soil nature, moisture and roughness. These layers are then analyzed in a GIS and supplemented if necessary on a scale of work that is based on different levels of operation. All these parameters were used to predict the vehicle performances, including speed, by comparing them to a database compiling standard cases of mobility. The resulting map shows the places where the vehicle can go without problems (green), with a relatively slow velocity (orange) and where it can't go (red) with a resolution of 30 m thanks to the DTED2 DEM. Calculations are performed using a Ground Vehicle Interaction (GVI) processing module, both off road but also on road, using a VMAP2i model, enriched by additional fields necessary to the ECORS estimators such as the maximum slope and sinuosity, presence of water or ice on the coating, presence of fords, bridges or tunnels.

RESULTS

In the present work, the ECORS system was tested on different test sites in France and Morocco. From a general point of view, the system delivers results in good agreement with observations (Land use classes, soil type, moisture, etc) measured on the field. A systematic study has been realized over different areas featured by different geographical contexts. The rate of success ranges between 97, 56 and 78%, depending on the considered area. Basically, most of errors were generated by two kinds of factors: 1-Some little escarpments (5 -10 m) were underestimated because of a lack of resolution of the DEM (30 m) that smoothed these topographic effects. In consequence, such barriers were wrongly considered as passable for the different vehicles; 2-Land use classification suffered also of some attribution errors between forest and bush. These confusions lead to pronounce NOGO on places considered as forest that were bush in reality. 3-Vehicle speed was also altered by some difficulties for estimating realistic soil moisture. Indeed, this quantity being estimated from climatic considerations is generally slightly different from what is observed in reality, and thus can produce variations in the final diagnosis. Despite such limitations, DAMs where successfully validated; they offer a valuable information for planning operational missions over badly known areas.

CONCLUSION AND FUTURE PLANS

In the present work, a decision-making system was designed, developed and tested in order to estimate vehicles mobility in both cross-country and roads modes. The ECORS system is based on processing, interpretation and data-fusion processes from remote sensing data, existing databases, including e-sources. Computing algorithms include slope computation, moisture index, land use and soil classification, road description, etc. Results are compiled in decision-making maps where two levels of information are available: NOGO and velocity maps, depending on data exhaustiveness. Validation processes have been carried out for verifying the coherence of these maps. Results show a good agreement between prediction and field observations, except for terrain with high complexity or non-exhaustive data. As perspectives, new field validation will be carried out in the next few months, particularly for testing the system in zones representative of operational areas. Figure .2: DAM for a 4x4 vehicle computed over an area located in the north-eastern part of Morocco; green, yellow, red colors refer to go, slow-go and no-go velocity classes, on and off-road respectively; blue and grey parts represent the floodplains and urban zones respectively.

Figure. 1 :

 1 Figure.1: general workflow of the ECORS system; terrain parameters are indicated in green, calculations in yellow and Decision Aid Maps (DAMs) in orange; GVI: Ground Vehicle Interaction processing.

ACKNOWLEDGEMENTS

The project was funded by DGA and has support from operational military French Forces. The ECORS project Team is composed by the following French organisms and industrials Brgm, Onera, Geo212, Sertit and Car&D.


M A Sbai
email: a.sbai@brgm.fr

PARALLEL PROCESSING OF THREE-DIMENSIONAL FIELD-SCALE REACTIVE TRANSPORT APPLICATIONS

In order to achieve field-scale reactive transport simulations in a reasonable CPU time, a parallel computational model is developed for RTAFF a fully three dimensional reactive transport simulator. The MPICH portable implementation of the message passing interface standard is used for communication among nodes in a Beowulf cluster of PC's. Parallel computations are carried out in different subdomains mapped to each processor in the Local Area Network during the advection and dispersion steps of the time-split procedure. The domain decomposition method is implemented using row-wise and column-wise overlapping subdomains simultaneously for the mixed hybrid finite element kernel. Computations of the reaction step are performed using a newly developed load balancing algorithm. This algorithm maps cells into processors based on global convergence history analysis from the previous time step. The efficiency of this approach is applied on a demonstrative field scale example showing that the load balancing algorithm efficiently localizes and redistributes reaction fronts cells among the processors. Furthermore, this application shows the scalability of the overall parallel algorithm and the achievement of several-orders-of-magnitude speedup among a range of processors. Modeling results provide new insights into detailed subsurface reactive-transport processes which could be taken simultaneously by scientists and engineers in either research or field-scale applications.

INTRODUCTION

The shape and scope of multi-component reactive transport models changed dramatically over the past decade. Advances in several branches of subsurface hydrology, aqueous and surface geochemistry has enabled an increased sophistication in modeling capabilities. Therefore, modern computer codes should offer prospective possibilities for modeling field-scale multi-component flow and transport coupled with complex and heterogeneous geochemical reactions. However, simply feeding many chemical reactions from a large geochemical database into large scale models is beyond the capabilities of currently available single processor workstations or PC's. Available parallel computing methodologies and ressources are therefore the only relevant means to achieve such simulations in a Date: March 13, 2006. 1 reasonable timeframe. When it seems that the path to a parallel computational model is a natural and straightforward approach, designing and efficiently implementing a scalable parallel reactive transport model is a non trivial task. Indeed, in the framework of an operator time-split technique (SNIA or SIA approaches [START_REF] Steefel | Approaches to modeling of reactive transport in porous media[END_REF]) a rough distribution of the geochemistry workload over the entire available processors will not yield an impressive parallel speedup. Pioneering work by [START_REF] Wheeler | Efficient parallel computation of spatially heterogeneous geochemical reactive transport[END_REF] showed examples of dynamic computational roughness typical to heterogeneous geochemical transport. Therefore, innovative load balancing algorithms have to be designed for the distributed geochemical calculations task. At the other hand, partitioning algorithms and preconditioning of high performance linear solvers for unstructured meshes are not obvious tasks. A number of available scientific toolkits for scalable parallel programming [START_REF] Balay | Efficient Management of Parallelism in Object Oriented Numerical Software Libraries, Modern Software Tools in Scientific Computing[END_REF] provide a management framework for the inter-communication and data partitioning details. However, such approach was not suited for the work presented in this paper, because our primary focus was on speeding up the computational geochemistry part. Indeed, previous experience on some selected practical problems showed that more than 90% simulation time was exclusively devoted to mixed equilibrium kinetics calculations. Parallel domain decomposition of the flow and transport solvers, although functional at the time of this writing, is implemented only on ijk structured grids. Adoption of parallel computing technology brings advantages and constraints on its own right. The most limiting being hardware availability for a user's community. Luckily, modern computer technology allows to setup clusters of PC's from cheap commodity hardware delivering massive computing power. The most popular distributed architecture falling in this category are Beowulf clusters. These are based on open-source software packaged with currently available distributions of Linux operating systems. Thus, by targetting such architectures the computer package will run unmodified for systems accessible to many modelers, as well as in highly priced shared memory workstations available in a number of sites at the national level. The primary objective of this work is to provide parallel numerical technology allowing us to perform reactive transport applications with several tens thousands cells (≈ 50,000), more than a dozen of master chemical species, and many thousands of time steps (≈ 10,000). This corresponds hopefully to a total CPU time reduction from a couple of days computing time to few hours. Next sections will describe the steps taken in our initial efforts on model development, evaluation, and testing on the basis of a three-dimensional model for deep acid gas (a mixture of H 2 S and CO 2 ) injection in a heterogeneous carbonated reservoir.

PARTAFF MODEL DEVELOPMENT

The computer program PARTAFF (Parallel Reactive Transport and Fluid Flow) simulates multi-component, reactive solute transport in three-dimensional saturated aquifers and single phase reservoirs. It can be viewed as a multi-processor extension of RTAFF simulator [START_REF] Sbai | A numerical model for miscible displacement of multicomponent reactive species[END_REF]. It is based on a modular structure including a simulation coordinator module, a fully implicit mixed hybrid finite element module for computing the flow in quadrilateral and hexahedral structured meshes, and an explicit high-resolution finite volume transport module. Meanwhile, geochemical calculations are based on a mixed formulation which couples a Gibbs free energy minimizer for constrained thermodynamic equilibrium equations with several integration methods for kinetically controlled reactions. Numerical solutions for pressure heads, Darcy's velocities, concentrations or activities of total dissolved components, mineral masses fractions, activities of user-selected secondary species, and rates of kinetically controlled reactions are obtained using a sequential approach. A four-stage advective-dispersive-reactive (ADR) operator-splitting algorithm [START_REF] Sbai | A numerical model for miscible displacement of multicomponent reactive species[END_REF] is used at each time step to decouple reaction and transport parts. Recently, a two-step time-split algorithm based on SIA [START_REF] Steefel | Approaches to modeling of reactive transport in porous media[END_REF] approach has been also incorporated in the simulator to minimize time-splitting errors. More details on the program formulation and validations should be found elsewhere [START_REF] Sbai | Numerical solution of mixed thermodynamic equilibrium and kinetically controlled equations in RTAFF code. I-evaluation of coupled geochemical transport[END_REF].

2.1. Hardware and software requirements. PARTAFF code use the portable implementation of the message passing interface standard provided by the MPICH project [START_REF] Gropp | A high-performance, portable implementation of the MPI message passing interface standard[END_REF]. We used only 29 routines from the 217 forming the whole library. It is a single program multiple data model (SPMD) that creates a fixed number of parallel tasks at the beginning of the program execution. Thus, a copy of the simulator is executed in each processor, while the data is being distributed over all available processors. This methodology does not allow for a smooth parallelization effort as provided by other standards (such as OpenMP [START_REF] Chandra | Parallel Programming in OpenMP[END_REF]) which may not pay-off for porting many existing models. This is not the case for PARTAFF which has been designed from the assets as a parallel and multiple platform code. Model evaluation and testing presented in this paper are performed on a Beowulf cluster of 12 dual-processor nodes with Intel Xeon processors running at 3.06 GHz and 2 GB local memory. All nodes use the Red Hat Linux operating system. The code was compiled with Portland Group high-performance suite of Fortran 90 and C/C++ compilers.

Parallel algorithms.

2.2.1. Domain decomposition (DD) of flow and transport computations. The simulation is carried out for 3D overlapping domains. The use of DD with overlapping [START_REF] Gropp | Domain decomposition methods in computational fluid dynamics[END_REF] layers have been recognized as being able to accelerate the convergence of the parallel computation. Hence, in this study the DD domain is decomposed into a number of overlapping subtasks, the computation of each subtask is carried out in one processor. In this study the domain is decomposed into a number of horizontal, transversal and vertical overlapping slabs using task decomposition as shown in Figure 1. The vertical decomposition is only allowed if no production or injection wells exist in the computational domain. The governing equations are solved for each subtask with appropriate boundary conditions. The solution process continues until the parallel preconditioned conjugate gradient iteration converge. Primary variables have to communicate between the processors at ghost points using blocking communication and the new values of primary variables are used in each subtask.

Load balancing of geochemical computations (The Heavy Light Hierarchically Randomized Splitting Algorithm).

A successful load balancing algorithm in this context should capture dynamic physics of the underlying problem. It is expected that the geochemistry workload will be highly dependent upon positions in space and time. As the advective and dispersive waves travel over the computational domain, zones of different behavior will develop throughout. Cells located in zones of nearly constant primary dependent variables [START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF] (total component concentrations) will need much less time to converge since the expected changes are minor. However, cells crossed by the traveling waves will exhibit the most significant workload. The steps involved in the developed algorithm for this task are given as follows:

(1) Estimate the required work for each grid-block assigned to a Beowulf node at the current time step. We assume that it equals the global geochemistry work calculated in the previous time step. A list of inactive cells undergoing minor total concentrations, C k , change such that:

max Nc k=1 | ∆C m-1 k |< ε 2 (1)
is updated in this step. All cells in this list do not contribute to the estimated work for processors, and are consequently excluded from the following treatments.

In Equation 1, m is the time stepping level, k is the chemical component index, N c is total number of components, ∆C m-1 k is the change of concentration of the k th component at the previous time level, and ε is the effective zero concentration (i.e. ε = 10 -30 M).

(2) Build a randomized spectral workload in Log 10 scale that would equilibrate the global work among the processors. This is done by constructing several subgroups of cells ranging from the most lighter to the most heavier, and assigning an average workload for each subgroup. Then, groups of cells are sent (through non blocking point to point MPICH calls) from heavy to light (listening) nodes according to the new average workload distribution.

(3) Perform parallel batch geochemical calculations for all nodes. First, local cells to each processor will be processed. Next, non local cells buffers moved in the previous step are processed. Finally, lighter nodes return newly computed buffers to sending processors to ensure proper updating with regard to the succeeding local transport step of the global time-split procedure. This algorithm does not only try to equilibrate the workload between the available processors. It also eliminates very light cells needing fewer iterations counts and satisfying equation 1 from the computational process. For typical applications, those noisy cells increase the parallel overhead, although not dramatically, and could be safely ignored letting the algorithm concentrating on space regions with noticeable concentrations changes. This improves the coarsening of the parallel algorithm, minimizes communication calls, and speeds up the cells classification process.

MODEL TESTING AND PERFORMANCE RESULTS

3.1. Description of the test case problem. The geochemistry behind this test case problem is taken from the study by [START_REF] Azaroual | Thermo-kinetic modelling of mass-exchange mechanisms between water-rock-gas induced by acid gas injection in carbonate reservoirs[END_REF] guided in the French Geological Survey. Its main purpose was to assess the chemical disequilibrium and sensitivities associated with heterogeneous kinetic rate laws typical for acid gas injection operations in a depleted carbonated reservoir. Batch and one-dimensional simulations were performed at the laboratory core scale (≈ 0.1m length). Herein a larger scale model is constructed to study the plume displacement mechanisms and the impact on rock mineralogy at the injector vicinity and farther away as well. The three dimensional model domain is a 1000m x 1000m x 100m box (101x101x5 uniform grid cells in the x, y and z directions respectively) where the mean vertical level is 1500m deep. A log-normally stochastic hydraulic conductivity field is constructed using the turning bands algorithm by [START_REF] Tompson | Implementation of the Three-Dimensional Turning Bands Random Field Generator[END_REF], with an average of 2x10 -2 m/day, dimensionless variance of 1.0, and correlation lengths equal to 5 times the grid cells size in the three principle dimensions. The injection flow rate is assumed to be continuous for the first 20 years of the total 50 years simulation period. The well is centered in the study area and is assumed to have a total penetration along the reservoir thickness.

Evaluation of model performance.

A methodological assessment of any parallel model performance involve a number of metrics. Common metrics are program relative speedup and relative efficiency defined respectively as [START_REF] Foster | Designing and Building Parallel Programs[END_REF]:

S n = T 1 T n (2)
E n (%) = 100 * S n n (3)
where T 1 and T n are run times with one and n processors, respectively. S n , the relative speedup with n processors, is a practical measure of the reduction in time achieved when increasing the number of processors. It is also as the number of apparent processors to the parallel computational engine. It is only limited by the number of uniprocessor chunks in the program and by the communication workload. E n , the relative efficiency with n processors, measures to which extent the parallel algorithm is effective during program execution.

Total execution time of the program is splitted into (i) CPU time in flow solver, (ii) CPU time in transport solver, (iii) CPU time in geochemical calculations, (iv) idle time, (v) communication time, and (vi) load balancing and parallel overhead time. Next, we will be interested only on relative speedup and relative efficiency associated to the transport and geochemical reaction tasks of the whole simulation. The load balancing algorithm takes only 0.3% of the total simulation time indicating that the geochemical load balancing algorithm is a coarse grain parallel task with a large workload to communication times ratio. Figure 2 illustrates the relative speedup and relative efficiency of the batch geochemical engine and solute transport models respectively. Strong testing of the problem scalability is performed as a function of increasing powers of 2 number of processors up to 16 processors limit. Clearly, it is shown that both parts of the simulation scale very well on the considered range of processors. The speedup is near to optimal up to eight processors, but maintains an important speedup of 13.45 for the 16 processor simulation. However, the transport speedup decreases more rapidly than the geochemical speedup. This is due to the increased communication overhead but also to the good scalability of the geochemistry workload making the transport speedup less apparent. However, this problem becomes of minor impact when much less CPU time is spent in the transport engine. For the range of processors taken into account the geochemical model efficiency remains always higher than 81%. The efficiency of the transport decreases much rapidly, even if it exceeds 70%. It is expected to be much less optimal for a high number of processors, deserving that much research must be guided for solving higher resolution problems on massively parallel computers in the future.

SUMMARY

This work presents a parallel simulation technology for multi-dimensional reactive transport applications. The technique involves a unified and portable set of parallel algorithm tools for saturated groundwater flow, solute transport, and geochemical calculations using MPICH message passing interface standard. This provides a simulation framework for systems ranging from multi-processor workstations, to cluster of PC's or workstations, and modern mainframe supercomputers. To demonstrate the relevance of this method, a test example for acid gas injection in a depleted carbonated reservoir is provided, and yields the following conclusions:

• The parallel algorithms are effective in reducing the total CPU time for field-scale reactive transport applications to an acceptable timeframe. • The relative importance of efficiently rebalancing the geochemistry workload is demonstrated, and shows to be critical to the scalability of the overall parallel computational model. • Presented parallel computing techniques can particularly enhance modeling capabilities, achieving several orders-of-magnitude speedup for field-scale and highresolution applications. • Much research must be conducted to evaluate the performance of the current algorithms in scenarios with increased geochemical complexity.

Figure 1 .

 1 Figure 1. Example of a 3D mesh partionned into eight (2x2x2) logical sub-domain processors. Ghost cells of each sub-domain (in grey shaded color) enable communications during preconditioned linear solver iterations.

Figure 2 .

 2 Figure 2. Relative speedup and relative efficiency of parallel transport and geochemical calculations for the test case problem.

ACKNOWLEDGMENTS

The author would like to thank Dr Mohamed Azaroual for many enlightening discussions on modeling geochemical transport processes.