Jocelyne Erhel
email: jocelyne.erhel@inria.fr

Souhila Sabit
email: sabitsouhila@yahoo.fr

Analysis of a global reactive transport model and results for the MoMaS benchmark

Keywords: reactive transport, global method, MoMaS benchmark

Reactive transport models are very useful for groundwater studies such as water quality, safety analysis of waste disposal, remediation, and so on. The MoMaS group dened a benchmark with several test cases. We present results obtained with a global method and show through these results the eciency of our numerical model.

Introduction

It is quite challenging to develop a numerical model for deep storage of nuclear waste. The time interval is very large (several thousands years), models are coupled and simulations must be accurate enough to be used for risk assessment. In most cases, chemistry must be included in models of deep geological storage. In addition to radioactive decay, chemical phenomna are numerous and include aqueous reactions, oxydo-reduction reactions, precipitation and dissolution reactions, ions exchanges, surface exchanges. These reactions can be either kinetic or at equilibrium. Models must handle species which are in groundwater systems and take into 10 account the mobile property of these species. It is thus necessary to consider a coupled model, where chemistry equations and radioactive decay are combined with transport of contaminants. These models are partial dierential equations (for transport, one equation for each species), and algebraic or dierential equations (for chemistry, a system at each grid point).

The MoMaS benchmark was designed as a set of academic examples, in order to run experiments with several methods and software [START_REF] Carrayrou | Looking for some reference solutions for the reactive transport benchmark of momas with specy[END_REF], [START_REF] Carrayrou | Reactive transport benchmark of momas[END_REF]. Several authors participated in the exercise [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF], [START_REF] De Dieuleveult | A global approach to reactive transport : application to the momas benchmark[END_REF], [START_REF] Homann | A parallel global-implicit 2-d solver for reactive transport problems in porous media based on a reduction scheme and its application to the momas benchmark problem[END_REF], [START_REF] Lagneau | Hytec results of the momas reactive transport benchmark[END_REF], [START_REF] Mayer | Solution of the momas reactive transport benchmark with min3p -model formulation and simulation results[END_REF]. A synthetic comparison of their results indicate that, for this benchmark, the fastest results were obtained with global approaches [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly non-linear and heterogeneous reactive transport problems. the momas benchmark case[END_REF].

In this paper, we show original results obtained with a global approach for the so-called 2D easy test case of the MoMaS benchmark. The model proposed in MoMaS is based on the introduction of total analytical concentrations, thanks to the linearity of the transport equation. It is a set of Partial Algebraic Dierential Equations. We use the method developed in [START_REF] De Dieuleveult | Un modèle numérique global et performant pour le couplage géochimie-transport[END_REF][START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF] and improved in [START_REF] Erhel | Computational science, engineering and technology series[END_REF][START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF], where we rst discretize in space, using a Finite Dierence scheme, then discretize in time, using an implicit multistep scheme, of the BDF family (Backward Dierentiation Formula). Thus at each time step, we have to solve nonlinear equations which we keep coupled.

Compared with [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF], three improvements were brought in [START_REF] Erhel | Computational science, engineering and technology series[END_REF][START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF][START_REF] Erhel | A global reactive transport model applied to the momas benchmark[END_REF]. First, a substitution technique, similar to global DSA methods, allows reducing the size of nonlinear or linear systems. The dierential variables are kept in the semi-discrete system, in order to use adaptive time steps and adaptive Jacobian updates. Logarithmic variables are very convenient to ensure the positivity of the concentrations and to compute the derivatives, but they can lead to severely ill-conditioned Jacobian matrices [START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF]. Thus, the second improvement consists in using non logarithmic variables, at the price of ensuring positivity during nonlinear iterations. In the benchmark, the rst component is inert, so the third improvement is to remove this component from the coupled equations.

The paper is organized as follows. In section 2, we describe the mathematical model of the MoMaS test case studied and provide our simulation results for ow and transport of the ve main species. The numerical method is dened in section 3, as well as three versions of our software GRT3D [START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF]. Finally, we analyse in section 4 the eciency of the improvements brought to the original global method. We discuss these results and outline future work in the concluding section 5.

MoMaS benchmark

The MoMaS group studies mathematical models and numerical simulations for nuclear waste disposal. A set of test cases were dened for transport reactive problems [START_REF] Carrayrou | Reactive transport benchmark of momas[END_REF]. Here, we make experiments with the so-called easy test case, in 2D. The computational domain is a rectangle with two porous media, see Figure 1. All dimensions are normalized, with length unit L and time unit T .

Flow simulations

The benchmark considers a steady saturated one-phase ow, with no source term. Flow is governed by Darcy's law and mass conservation, giving mathe-matical equations [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF] where the pressure h and the Darcy's velocity q are the unknowns and K is the hydraulic conductivity.

   q = -K∇h, ∇q = 0. (1)
Darcy's velocity is related to the porosity ε and to the pore velocity v by q = εv.

Medium A has a high conductivity and a low porosity, whereas medium B has a low conductivity and a high porosity, see Table 1.

Medium A Medium B Porosity ε 0.25 0.5 Conductivity K (L.T -1 ) 10 -2 10 -5
Table 1: Flow conditions for MoMaS benchmark.

Boundary conditions must be prescribed to complete the PDE system. At outow, the pressure is given by h = 1L. At both inows, the ow velocity is given by q = 2.25 × 10 -2 L.T -1 , whereas no ow condition is applied at other boundaries.

We simulate the ow equations with the software MODFLOW [START_REF] Harbaugh | MODFLOW-2000, the U.S. Geological Survey modular ground-water model User guide to modularization concepts and the Ground-Water Flow Process[END_REF]. Figures 2 et 3 show the velocity and the pressure computed with a mesh of 40 × 84 cells [START_REF] De Dieuleveult | Un modèle numérique global et performant pour le couplage géochimie-transport[END_REF], [START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF].

Transport simulations

Transport of mass in groundwater is governed by advection and dispersion.

There is no source term in the benchmark test. The concentration c of an inert solute is the solution of the PDE (2) expressing a mass conservation law.

ε ∂c ∂t = ∇.(D ∇c) -∇.(q c), (2)
where the dispersion tensor D is given by We consider the advective test case of the MoMaS benchmark, without molecular diusion, see Table 2.

D = εd m I + α T q I + (α L -α T ) qq T q .
Initial conditions at time t = 0T are applied to the geochemistry system, see The rst chemical component is a spectator ion, which behaves like a inert solute. For this component, we can simulate the transport equations with the software MT3D [START_REF] Zheng | MT3DMS : a modular three-dimensional multi-species model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems : documentation and user's guide[END_REF]. Figure 4 represents the concentration at various times, obtained with a ne mesh of 80 × 168 cells.

Geochemistry simulations

In the MoMaS test case, the geochemical system has N c = 4 aqueous components c j and N s = 1 xed component s j . They react with N α = 5 aqueous secondary species α i and N β = 2 xed secondary species β i . There is no precipitation dissolution. All the coecients of activities are equal to 1 and the variables c j , s j , α i , β i are the concentrations of the species.

Stoichiometric coecients and equilibrium constants are given in Table 3. In the original benchmark, the constant K α5 is equal to 10 +35 but we replace it by 10 +6 . Indeed, with such a very large constant, the chemical nonlinear system is highly ill-conditioned and inaccuracies may appear for some values of c and s.

For example, we succeeded to run 1D simulations with K α5 = 10 +35 [18], but we experienced numerical diculties in 2D.

Table 3 is summarized in Table 4 by using algebraic notations, with the matrices S, A, B and the vectors c, s, α, β, K α , K s .

Secondary species can be computed thanks to the mass action laws (3) describing the chemical reactions. Moreover, the concentrations of the components

c 1 c 2 c 3 c 4 s K α 1 0 -1 0 0 0 10 -12 α 2 0 1 1 0 0 1 α 3 0 -1 0 1 0 1 α 4 0 -4 1 3 0 0.1 α 5 0 4 3 1 0 10 +6 β 1 0 3 1 0 1 10 +6 β 2 0 -3 0 1 2 10 -1
Table 3: Stoichiometric coecients and equilibrium constants for MoMaS benchmark (with

K α5 modied). c s K α S 0 K α β A B K s
Table 4: Algebraic representation of Table 3. must be positive.

               α i (c) = K αi Nc j=1 c Sij j , i = 1, . . . N α , β i (c, s) = K βi Nc j=1 c Aij j Ns j=1 s Bij j , i = 1, . . . N β , c j ≥ 0, j = 1, . . . N c , s j ≥ 0, j = 1, . . . N s . (3)
Then the model applies the mass conservation law to the chemical system.

   c + S T α(c) + A T β(c, s) = T, s + B T β(c, s) = W, (4)
where T and W are respectively the total analytical concentrations for mobile and xed components. In a closed system, these quantities are known, but they vary in time and space when the species are transported by water.

Reactive transport

Transport is now governed not only by advection and dispersion, but also by reaction. A mass conservation equation can be written for each component and each secondary species, where the reaction terms are unknowns in the context of reactions at equilibrium. However, because the dispersion tensor is the same for all the species and the transport terms are linear, it is possible to compute a linear combination of all the equations. This ends up to mass conservation equations ( 5) applied to the total analytical concentrations T , where the reaction terms cancel and disappear. Moreover, the total analytical concentrations W remain constant and given.

ε ∂T j ∂t = ∇.(D ∇C j (c)) -∇.(q C j (c)), j = 1, . . . N c , (5)
where C(c) = c + S T α(c) is the total mobile concentration.

Equations ( 5) are coupled with the mass action laws (3) and the mass conservation laws (4) written at each point of the computational domain.

Initial conditions for the 5 components are given in Table 5. It can be noted that T 3 = 0 implies that c 3 = α 2 = α 4 = α 5 = β 1 = 0, because all the associated stoichiometric coecients are strictly positive. Therefore, the equilibrium constants K α2 , K α4 , K α5 , K β1 have no mathematical eect on the initial equilibrium. Thus the original initial conditions of components remained the same when we changed the value of K α5 .

Boundary conditions at both inows are given in Table 6 for the injection and leaching periods. Again, it can be noted that T 3 = 0 during the leaching period, and T 4 = 0 during the injection period. Since the stoichiometric coecients associated to c 4 are also strictly positive, c 4 = α 3 = α 4 = α 5 = β 2 = 0 during the injection period. Thus, at inow, the equilibrium constants K α4 , K α5 , K β2 have no mathematical eect on the equilibrium.

Figure 5 shows the concentrations of the aqueous components c 4 at various times, using a mesh of size 80 × 168. They cannot be compared exactly to the results of the literature [START_REF] Carrayrou | Looking for some reference solutions for the reactive transport benchmark of momas with specy[END_REF][START_REF] Homann | A parallel global-implicit 2-d solver for reactive transport problems in porous media based on a reduction scheme and its application to the momas benchmark problem[END_REF][START_REF] Lagneau | Hytec results of the momas reactive transport benchmark[END_REF][START_REF] Mayer | Solution of the momas reactive transport benchmark with min3p -model formulation and simulation results[END_REF], because we changed the constant K α5 , but they are very close. More results can be found in [START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF].

T 1 T 2 T 3 T 4 W Medium A 0 -2 0 2 1 Medium B 0 -2 0 2 10

Numerical method

Equations ( 3),( 4),( 5) form a PDAE system composed of algebraic and PDE equations. We assume that it is well-posed.

Following a method of lines, we rst discretize in space, using a nite difference scheme implemented in the software MT3D. We dene a regular mesh of the domain with a rectangle of N x × N y cells. We get a DAE semi-discrete system, where the algebraic chemistry equations are written at each point of the mesh [START_REF] De Dieuleveult | Un modèle numérique global et performant pour le couplage géochimie-transport[END_REF][START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF].

Because the system is sti, we use an implicit scheme, involving a global nonlinear system at each time step, coupling transport and chemistry equations.

We use a BDF scheme implemented in the software SUNDIALS [START_REF] Hindmarsh | SUNDIALS : Suite of nonlinear and dierential/algebraic equation solvers[END_REF]. At each time step, the nonlinear system is solved with a Newton's method. We provide the function of the DAE equation and its derivative, whereas the module IDA of SUNDIALS provides the implicit scheme.

Time discretization comes with an adaptive time step, which controls both the accuracy of the approximation and the convergence of Newton's iterations.

It allows choosing large time steps when possible and saves CPU time.

We solve the linearized equations with a sparse direct solver, implemented in the software library UMFPACK [START_REF] Davis | Algorithm 832 : Umfpack, an unsymmetric-pattern multifrontal method[END_REF], which we have interfaced with SUN-DIALS. We thus factorize the Jacobian matrix and use the triangular factors to solve the linearized equations. The software keeps the Jacobian of linearized Newton's iterations frozen while convergence is fast enough, saving updates and factorizations thus CPU time.

We implemented three versions of our method in the software suite GRT3D [START_REF] Sabit | Suite Logicielle GRT3D (Global Reactive Transport 3D), Rapport de contrat Andra[END_REF][START_REF] Erhel | Computational science, engineering and technology series[END_REF][START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF]. Mass action laws can be linearized by introducing logarithmic variables, if they are strictly positive. We introduce ĉ = log(c), ŝ = log(s), so that equation ( 3) becomes (written with matrices and vectors) :

   α(ĉ) = exp (log K α + S ĉ) , β(ĉ, ŝ) = exp (log K β + A ĉ + B ŝ) ,
and c, s are replaced by exp(ĉ), exp(ŝ) in equation ( 4).

In a rst version, called GRT3D, we used as primary variables the total T , the total mobile C and the components concentrations ĉ, ŝ, so that we had (3N c + N s ) variables at each cell. In a second version, called GRT3DRL, we reduced the size of linearized systems by using a substitution approach in the linear systems to be solved. We eliminated the variables T and C and kept only the variables ĉ, ŝ, thus (N c +N s ) variables at each cell. This technique saves CPU time by factorizing a smaller matrix. In a third version, called GRT3DRSL, we used the same algorithms, but with the non logarithmic variables c, s. We also eliminated the variables T and C in the linearized equations. In this case, we have to ensure the positivity of concentrations. This is done in Newton's method, where approximate concentrations are forced to be positive at each iteration.

Performance analysis

In this section, we analyze the accuracy and the computational time of our simulations, using the three versions of our software. Initial and boundary condi- Here we show a selection of the results which can be found in [START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF]. For each component, we compare the concentrations obtained with the three versions of the software, keeping or removing the inert component. The algorithms are mathematically equivalent but not numerically, therefore we compare their accuracy. We take as reference the result obtained with GRT3DRSL, that is the version with substitution and without logarithms. For each component c i or s i , we compute the dierence between this reference and the results obtained with other versions, in the following way :

E ci = 1 N m N t n,j (c i (m j , t n ) -c i (m j , t n )) 2 , ( 6
)
where m j is the cell number and t n is the timestep number, with N m cells and N t timesteps, ci is the reference concentration and c i is the other concentration.

Table 7 shows the deviations for the component c 2 , those for the other components can be found in [START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF]. Since these dierences are quite small for the three meshes and for all versions, we conclude that we do not introduce numerical artefacts by using substitution or by using logarithms. Table 8: System size using GRT3D and GRT3DRL with and without the rst inert component.

The reduction of the system size is given in Table 8. If we keep the inert component and apply substitution, the system size is reduced from Table 9: CPU time of the three versions with the rst inert component.

(3N c + N s ) × N m = 13N m in GRT3D to (N c + N s ) × N m = 5N m in
Clearly, reducing the system size has a direct impact on the CPU time. The total CPU time for the three versions is given in Table 9 for computations with the rst inert component and in CPU time decreases also when variables are not logarithmic. In order to analyze these eects, we use the measures taken at each external timestep of the software SUNDIALS. In Figure 6, we plot the CPU time versus the simulated time, for the six simulations (three versions of software, with and without the rst component), and for the ne mesh. During all the simulated time interval, the substitution approach, as well as the elimination of the rst component, reduce the CPU time. Also, the use of non logarithmic variables is much more ecient.

For all the simulations, the CPU time increases rapidly until about t = 1000T then it levels o until time t = 5000T and increases again rapidly when the inow boundary conditions change. We analyse the outputs of the software IDA in order to explain this behavior.

Eect of adaptive time step

The IDA solver in the software SUNDIALS adjusts automatically the time step, in order to control both the accuracy of time discretization and the convergence of Newton's iterations. In Figure 7, we plot the number of time steps taken during the whole simulation. Here, we show only results without the rst component and for the ne mesh, since other results are very similar.

As expected, the CPU time is strongly correlated with the number of time steps. The initial and boundary conditions are dicult to handle and require small time steps during a period of time until larger time steps can be taken.

These results demonstrate the eciency of adaptive time step.

Eect of modied Newton's iterations

The number of time steps has a direct impact on the number of Newton's iterations which can be measured by the number of linear system solvings. It is also correlated to the number of updates of the Jacobian matrix, which can be measured by the number of matrix factorizations. Indeed, in our simulations, we use a direct sparse linear solver implemented in the software UMFPACK.

Thus each time the Jacobian is updated, it has to be factorized in order to solve the linear systems associated to this new Jacobian.

In Figure 8, we plot the number of linear solvings and matrix factorizations, for the three versions and for the ne mesh. Theses numbers are in logarithmic scale. Here, the rst component is not considered, similar results are obtained when including it.

Again, we observe a sharp increase of these numbers at the beginning of the simulation and a small jump after leaching. At time t = 0, the number of solvings as well as the number of factorizations are slightly smaller when reducing the system size, from GRT3D to GRT3DRL. They are clearly smaller when using non logarithmic variables. This result is probably due to a better condition number of the Jacobian matrix at time t = 0.

For the three versions, the number of matrix factorizations is much smaller than the number of linear system solvings. This result demonstrates the eciency of adaptive modied Newton's iterations. Indeed, the Jacobian is kept frozen during several time steps and during nonlinear iterations, so that factorization does not occur as often as solving.

Algorithmic complexity of Newton's iterations

It is well-known that the algorithmic complexity of sparse matrix factorization is much higher than the complexity of sparse system solving, which involves sparse triangular matrices. Thus, in order to fully measure the eciency of the adaptive update of the Jacobian matrix, we have to analyze the CPU time of the Newton's iterations.

In Figure 9, we plot the CPU time (in logarithmic scale) of the linear solvings and matrix factorizations, for the three versions and the ne mesh. We see again the eect of the numerical diculties at time t = 0 and the dierences between the three versions.

Whereas the number of factorizations is smaller than the number of solvings, the CPU time does not behave in the same way. For the small and intermediate meshes, factorization time is smaller than solving time but it becomes larger than solving time for the ne mesh. Because the algorithmic complexity of the factorization is higher, the CPU time increases much faster when the mesh is rened.

This is also illustrated by Figure 10, where we plot the percentage of time spent in linear steps in function of the mesh size, measured by the number of cells. These percentages are computed versus the total CPU time, at the end of the simulations. Computations are done here with the version without logarithms (GRT3DRSL) and without the rst component c 1 .

Factorization time becomes relatively more important as the mesh size increases and eventually the lines intersect roughly at the size of the ne mesh. It can also be noted that the addition of both percentages increases with the mesh size. For the ne mesh, more than 90% of CPU time is spent in the Newton's iterations. Therefore, it is really important to reduce the number of iterations and their computational cost.

Concluding remarks

In this paper, we have proposed a global approach for simulating reactive transport equations, where transport is coupled with geochemistry. The global method is based on a DAE formulation of the semi-discrete system obtained after space discretization. An implicit time scheme ensures stability and provides an adaptive time step with an adaptive update of the Jacobian matrix. This feature is very ecient as illustrated in our numerical results with the MoMaS benchmark. In order to use this adaptive time scheme, we keep the dierential variables, but we eliminate them at the linear level, thus reducing the system size. This is also very ecient as demonstrated in our results. We could also apply the substitution at the nonlinear level [START_REF] Erhel | Computational science, engineering and technology series[END_REF][START_REF] Sabit | Les méthodes numériques de transport réactif[END_REF]. In the MoMaS benchmark, some concentrations are set to 0 at the intial time and at inow boundaries.

When using logarithmic variables and very small initial concentrations instead of 0, this leads to ill-conditioned systems. It appears that it is more ecient to use nonlogarithmic variables for the test case studied here.

We use a sparse direct solver which rst factorizes the Jacobian matrix, then solves two triangular systems. This second step has a much smaller complexity than the rst one. Thanks to the adaptive update of the Jacobian, the computational cost of the factorization remains low for the meshes considered here.

However, it is clear that this cost will eventually dominate for larger systems.

Therefore, we investigate parallel iterative solvers in order to tackle 3D problems [START_REF] Pacull | Parallel preconditioners for 3d global reactive transport[END_REF].

The test case considered here denes chemistry at equilibrium, without precipitation or dissolution. These reactions are very challenging because minerals can appear and disappear. We study some mathematical issues of this problem [START_REF] Migot | Analyse mathématique de modèles géochimiques[END_REF]. We also plan to include kinetic reactions.

Figure 1 :

 1 Figure 1: Computational domain for MoMaS Benchmark (dimensions are in the unit length L).

Figure 2 :

 2 Figure 2: Velocity eld with a mesh of 40 × 84 cells.

Figure 3 :

 3 Figure 3: Pressure head with a mesh of 40 × 84 cells.

tions are changed from 0 to 10 -

 10 10 when necessary for logarithmic variables (GRT3D and GRT3DRL). All experiments are done in sequential, on a Intel Xeon computer, with 24 MB of RAM and with 12 MB of cache memory. The rectangle is meshed by a regular grid of N m = N x × N y cells. Simulations are done with three meshes : a coarse mesh with N m = 20 × 42 cells, an intermediate mesh with N m = 40×84 cells, and a ne mesh with N m = 80×168 cells.

  GRT3DRL and GRT3DRSL. If we remove the inert component, we get (N c -1) aqueous components, and obtain a reduction thanks to substitution from 10N m to 4N m .

  (a) Time t = 1T (b) Time t = 4000T (c) Time t = 5010T (d) Time t = 5800T

Figure 4 :

 4 Figure 4: Concentrations of the inert solute c 1 at dierent times (with dierent scales).Simulations are done with a mesh of size 80 × 168.

20 (

 20 a) Time t = 1T (b) Time t = 4000T (c) Time t = 5010T (d) Time t = 5800T

Figure 5 :

 5 Figure 5: Concentrations of the solute c 4 at dierent times (with dierent scales). Simulations are done with a mesh of size 80 × 168.

Figure 6 :

 6 Figure 6: CPU time during the simulation for the six variants.

Figure 7 :

 7 Figure 7: Number of time steps during the simulation.

Figure 8 :

 8 Figure 8: Number of linear solvings and matrix factorizations during the simulation (logarithmic scale).

Figure 9 :

 9 Figure 9: CPU time during the simulation of the linear system solvings and the matrix factorizations (logarithmic scale).

Figure 10 :

 10 Figure 10: Percentage of CPU time in linear solvings and factorizations versus the mesh size using GRT3DRSL.

Table 5 .

 5 The nal time is t = 6000T .

	At impermeable boundaries, a no total ux condition is imposed. At ouow,
	a zero concentration gradient is imposed. At both inows, concentration is pres-
	cribed, with values given in Table 6. Injection occurs during a rst period of
	time until t = 5000T , followed by a leaching period until the end.

Table 2 :

 2 Transport conditions for MoMaS benchmark.

Table 5 :

 5 Initial conditions for MoMaS benchmark.T 1 T 2 T 3 T 4 W

	Injection t ∈ [0, 5000]	0.3 0.3 0.3 0	0
	Leaching t ∈ [5000, 6000]	0	-2	0	2	0

Table 6 :

 6 Inow conditions for MoMaS benchmark.

Table 7 :

 7 Comparison of accuracy for the component c

	4.1. Accuracy
	Mesh with the rst component without the rst component
	GRT3D	GRT3DRL	GRT3D	GRT3DRL
	20x42 6.402528E-08 3.112242E-09 3.227338E-08	2.568312E-09
	40x84 6.432335E-08 6.467349E-09 3.040179E-08	3.034366E-09
	80x168 5.802496E-08 5.286914E-09 6.091884E-08	5.727954E-09

[START_REF] Carrayrou | Looking for some reference solutions for the reactive transport benchmark of momas with specy[END_REF]

. The dierence with GRT3DRSL is computed from

[START_REF] De Dieuleveult | Un modèle numérique global et performant pour le couplage géochimie-transport[END_REF]

.

Table 10 :

 10 Table 10 without it. CPU time of the three versions without the rst inert component.

	Mesh	GRT3D	GRT3DRL	GRT3DRSL
	20x42	40 m 20 s	19 m 1 s	13 m 49 s
	40x84	7 h 20 m 22 s	4 h 4 m 1 s	2 h 12 m 45 s
	80x168 2 j 18 h 22 m 46 s 2 j 6 h 56 m 42 s 1 j 8 h 20 m 46 s

Acknowledgements

This work was partly funded by a grant from ANDRA and by a grant from ANR (H2MNO4 project).


Romain Pinquié
email: romain.pinquie@ensam.eu

Philippe Véron
email: philippe.veron@ensam.eu

Frédéric Segonds
email: frederic.segonds@ensam.eu

Nicolas Croué
email: nicolas.croue@keonys.com

Natural Language Processing of Requirements for Model-Based Product Design with ENOVIA/CATIA V6

Keywords: Functional digital mock-up, ENOVIA V6, CATIA V6, Natural language processing, Requirements, Parametric modelling

The enterprise level software application that supports the strategic product-centric, lifecycle-oriented and information-driven Product Lifecycle Management business approach should enable engineers to develop and manage requirements within a Functional Digital Mock-Up. The integrated, model-based product design ENOVIA/CATIA V6 RFLP environment makes it possible to use parametric modelling among requirements, functions, logical units and physical organs. Simulation can therefore be used to verify that the design artefacts comply with the requirements. Nevertheless, when dealing with document-based specifications, the definition of the knowledge parameters for each requirement is a labour-intensive task. Indeed, analysts have no other alternative than to go through the voluminous specifications to identify the values of the performance requirements and design constraints, and to translate them into knowledge parameters. We propose to use natural language processing techniques to automatically generate Parametric Property-Based Requirements from unstructured and semi-structured specifications. We illustrate our approach through the design of a mechanical ring.

In 1990, Gero [START_REF] Gero | Design prototypes: a knowledge representation schema for design[END_REF] proposed the FBS ontology where F stands for the set of functions, B for the set of expected behaviours (Be) and the set of actual behaviours (Bs), and S for the structure. In [START_REF] Christophe | RFBS: a model for knowledge representation of conceptual design[END_REF], Christophe extends the FBS ontology to RFBS by including the R for requirements. Back in the nineties, in his theory of axiomatic design, Suh [START_REF] Suh | Axiomatic design: advances and applications[END_REF] defined four domains of activities: the customer domain, the functional domain, the physical domain and the process domain. Stepping back and looking at these product design methods, which could also be assimilated to the systems engineering process [START_REF]Systems and software engineering -System life cycle processes[END_REF], we notice that product design relies upon an iterative process among requirements, functions, behaviours and structures. The Dassault Système's ENOVIA/CATIA V6 software solution proposes a similar integrated product design model named RFLP [START_REF] Kleiner | Model based design with systems engineering based on RFLP using V6[END_REF]. The R is for ENOVIA V6 Requirements, a requirements management workbench. The F, L and P layers are used to recursively break down the complexity of the design problem according to the Functional, Logical and Physical viewpoints of the product. This design approach follows from Descartes' reductionism method that consists in understanding a complicated problem by investigating simple parts and then reassembling each part to recreate the whole. In RFLP, the functional layer (F) relies upon a Functional Flow Block Diagram to design functional architectures in which functions transform material, energy or information input flows into output flows whose consistency is ensured by the matching of input and output typed-ports. The logical layer (L) is the behavioural viewpoint of the product and is materialised by a logical architecture within which each logical unit's behaviour is equation-based modelled with the Modelica1 language. Modelica models are executable thanks to the Dynamic Behaviour Modelling workbench that is the integration of Dymola2 within CATIA V6. Finally, the physical layer (P) is very similar to the CATIA V53 CAD modeller.

The integrated RFLP product design environment enables designers to define implementation links between a pair of requirements, functions, logical units or physical organs so as to trace implementation relationships thanks to a traceability matrix. In addition to the traceability capability, the tight integration of ENOVIA V6 Requirements and CATIA V6 offers parametric modelling functionalities that can be used to make sure that the design artefacts comply with the requirements.

Problematic

Among the product life cycle phases defined by Terzi [START_REF] Terzi | Product Lifecycle Managementfrom its history to its new role[END_REF], we focus on the requirement analysis phase without addressing the management of the requirements during the downstream detailed design and testing life cycle phases.

Nowadays, a set of requirements is usually very large. Indeed, with the ever-increasing complexity of products and their relentless customisation, the mushrooming accumulation of legal documents, let alone the geographically dispersed teams through whom products are developed, a supplier is faced with a staggering increase in the number of requirements. For instance, at Mercedes-Benz, the size of a system-of-interest (SOI) specification varies from 60 to 2000 pages and prescribes between 1000 and 50 000 requirements [START_REF] Houdek | Challenges in automotive requirements engineering[END_REF]. In addition to the massive volume of requirements, most specifications are unstructured documentse.g. Word, PDFand 79 % of requirements are written in unrestricted natural language [START_REF] Mich | Market research for requirements analysis using linguistic tools[END_REF].

For all these reasons, in a "buy approach" of a "make vs buy" decision, OEMs struggle to deliver products that comply with the legal and contractor's requirements. Indeed, when an OEM collects the specifications and the applicable documents the specification refers to, he has no other alternative than to go through the documents to identify the applicable requirements so as to, in fine, provide a product that complies with the contractor's requirements. There are four standard verification methods: inspection, analysis/simulation, demonstration, and test [START_REF]Systems and software engineering -Life cycle processes requirements engineering[END_REF]. In this paper, we benefit from the simulation method that the parametric modelling CATIA V6's capabilities offer in its integrated RFLP product design environment. Parametric modelling-based verification can be very time consuming since designers have to: (1) read the specifications, (2) identify the values of the performance requirements and the design constraints, (3) model the values of the performance requirements and the design constraints as requirements' knowledge parameters, (4) design the behavioural and structural artefacts using parametric modelling, and (5) define the knowledge verification rules that map requirements' knowledge parameters with design artefacts' knowledge parameters so as to verify their compliance.

In a "make approach" there is no exchange of document-based specification. Therefore, before designing, the company simultaneously prescribes the product's requirements and the requirements' knowledge parameters into ENOVIA V6 Requirements. However, in a "buy approach", the OEM has to move the requirements from the unstructured specification documents to the ENOVIA V6 Requirement database and manually build requirements' knowledge parameters.

Literature Review and Proposition

According to Lash [START_REF] Lash | Computational representation of linguistics semantics for requirements analysis in engineering design[END_REF], modal verbs such as shall, must and should are key lexical features for classifying sentences corresponding to requirement statements. In [START_REF] Coatanéa | Requirements models for engineering, procurement and interoperability: a graph and power laws vision of requirements engineering[END_REF], Coatanéa et al. use the Stanford Parser to apply natural language processing (NLP) techniques such as sentence splitting, tokenization and POS-tagging so as to create a binary rules-based classifier based on the presence or absence of a modal verb in a sentence. Zeni et al. [START_REF] Zeni | GaiusT: supporting the extraction of rights and obligations for regulatory compliance[END_REF] propose GaiuST, a framework that extracts legal requirements for ensuring regulatory compliance.

Few research studies attempt to extract text-based requirements (TBRs) from unstructured specifications; however, none of them address the challenge of extracting the values of performance requirements and design constraints to ease simulation-based design verification.

To avoid the very time-consuming requirements' knowledge parameters definition process, we propose a NLP pipeline to extract TBRs from unstructured and semi-structured specifications and to model them as Property-Based Requirements (PBRs). PBRs are used to automatically generate Parametric PBRs (PPRBs) in ENOVIA V6 Requirements. Finally, while designing with CATIA V6 RFLP, designers define behavioural and structural design knowledge parameters that are manually mapped to PPBRs thanks to parametric knowledge verification rules.

The model-based product design process that we present is twofold: (1) we extract TBRs from document-based specifications and transform them into PPBRs in ENOVIA V6 R2015X; (2) we exploit the PPBRs in an integrated, parametric, model-based product design synthesis with CATIA V6 R2015X.

From Unstructured Specifications to PPBRs

Before presenting the NLP pipeline that generates the PBRs, we must present the concepts of PBR and PPBR.

As Micouin introduces in [START_REF] Micouin | Toward a property based requirement theory: system requirements structured as a semilattice[END_REF], a PBR is an unambiguous formal definition of a requirement as a predicate and is defined as follows:

PBR : When C ! val O:P ð Þ 2 D ð1Þ
This formal statement means: "When the condition C is true, the property P of object type O is actual and its value shall belong to the domain D". A relevant characteristic of the concept of PBR is that it is grammar-free, i.e. a PBR does not have any particular syntactic structure and can therefore be implemented with various modelling language such as VHDL-AMS [START_REF] Micouin | Property-model methodology: a model-based systems engineering approach using VHDL-AMS[END_REF] and Modelica.

By combining the PBR theory with parametric CAD modelling, we coin the concept of Parametric PBR (PPBR). A PPBR is a PBR that is implemented with a parametric CAD modeller thanks to knowledge parameters and knowledge verification rules. In ENOVIA V6, a PPBR is analogous to the formal combination of an Object (O)the subject in the requirement statement attributewith one or several knowledge parameters that define the boundaries of the constrained domain (D) of a property (P), whereas the condition (C) is a CATIA V6 knowledge verification rule that is manually defined by designers while designing behavioural and structural artefacts.

The generation of PPBRs from TBRs requires: (1) an NLP 4 pipeline to generate an XML specification that stores a set of PBRs derived from TBRs, and (2) to interpret the XML specification of PBRs for creating PPBRs in ENOVIA V6 Requirements.

To derive PBRs from TBRs we implemented the following NLP pipeline:

Step 1 <Uploading>: The user uploads one or several specifications whose extension is .doc(x) (Word), .odf (OpenOffice), .pdf (Portable Document Format), .xls (x) (Excel), .xmi (SysML requirements diagram). While uploading, the file uploader item gets the input stream of each specification.

Step 2 <Parsing>: We trigger a specific parser according to the file extension of each specification. If it is a .doc or .odf, the parser uses the Apache Tika 5 API to extract each specification content and transform it into .html semi-structured data. We transform the content into HTML because it makes the analysis of tables, lists, headings, etc. a lot easier. The headings of .doc and .odf help us to identify the sections. Sections are used for multi-threading to run processing tasks in parallel. If the extension of the specification corresponds to a .pdf, we use the native capability of Word to convert from .pdf into .doc. Then, as for .doc, we use the Apache Tika API to convert from .doc into .html. Once we get the .html specification, we verify whether the .pdf was generated with Word or OpenOffice by looking for the header, footer or table HTML tags, which are not present in a .pdf that was created with another text editor such as LaTeX. If we find that the .pdf was generated with Word or OpenOffice, then we call the .doc parser; otherwise, we use the .pdf parser that relies upon the Apache Tika API and various regular expressions. The second scenario is less accurate as we lose the structures (tables, enumerations, footer, header, etc.). The .xls(x) parser uses the Apache POI 6 API to parse the textual content of each cell. We make the hypothesis that each cell is a sentence. Finally, the .xml parser extracts the requirement statements between the XML tags of a SysML requirement diagram.

Step 3 <Tokenization>: The Stanford CoreNLP [START_REF] Manning | The stanford CoreNLP natural language processing toolkit[END_REF] Tokenizer7 API iteratively tokenizes each specification content, that is, it chops the textual content up into pieces of a sequence of characters that are grouped together as a useful semantic unit for processing, the tokens. We store the tokens in a term-sentence matrix whose rows are sentences and columns are tokens that make up each sentence.

Step 4 <Lemmatization>: The Stanford CoreNLP Lemmatizer API iteratively normalises each token by removing the inflectional ending and returns the dictionary form, the lemma. For instance, lemmatization reduces the tokens "requires", "required" and "require" to their canonical form "require". This enables us to increase the recall in step 8 <Classification>.

Step 5 <POS-tagging>: The Stanford CoreNLP POS-tagger8 API iteratively POS tags each token, that is, it annotates each token with its grammatical category (noun, verb, adjective, adverb, etc.), the Part Of Speech (POS).

Step 6 <Sentence splitting>: The Stanford CoreNLP API iteratively splits each textual specification content into sentences.

Step 7 <Sentences cleaning>: We use various regular expressions and analyse HTML tags to clean the sentences. For instance, we rebuild sentences from enumerations, get rid of the headings, headers, footers and informative sections (introduction, scope, table of content, glossary, list of acronyms, etc.) that may generate false positives, and extract the content of .pdf tables.

Step 8 <Classification>: We use a knowledge engineeringa.k.a rules-basedtext classification approach [START_REF] Feldman | The Text Mining Handbook. Advanced Approaches in Analyzing Unstructured Data[END_REF] to binary classify each sentence into a "requirement" vs "non-requirement" class. The algorithm iteratively traverses the matrix whose rows are sentences and columns are lemmas. For each iteration, if the condition "token i of sentence j = a prescriptive term 2 {shall, must, should, have to, require, need, want, expect, wish or desire}" is true, the current sentence i is classified as a requirement.

Step 9 <Dependencies analysis>: The Stanford CoreNLP Dependencies Analyzer9 [START_REF] Cer | Parsing to Stanford dependencies: trade-offs between speed and accuracy[END_REF] API iteratively analyses each requirement to generate a semantic graph within which we identify the numeric dependencies and extract the source and target nodes of each dependency. The source of a numerical dependency is a numerical token annotated with the POS tag (CD), whereas the target is a word.

Step 10 <Classification>: While going through the dependencies list of each requirement, we check whether the word stored in the target node of each dependency is a physical unit such as N, °C, kg, Pa, etc. using a resource file that collects all existing physical units under its abbreviated and expanded forme.g. N and Newton. Each time a given numerical dependency is classified as a physical numerical dependency, we add a third attribute from our resource file that is the dimension of the physical unite.g. Force for the unit N or Newton.

Step 11 <PBR Pattern Analysis>: A well-written TBR prescribing a functional level of performance or a design constraint usually follows three distinct syntactic patterns (Pattern 1, 2 and 3) [START_REF]INCOSE: Guide for writing requirements. Requirements working group[END_REF]. Note that the condition is not always specified; consequently, there is one more syntactic pattern (Pattern 4).

Pattern 1: <Prescriptive> <Domain> <Condition> -PDC

The Control_Subsystem shall open the Inlet_Valve in less than 3 seconds when the temperature of water in the Boiler is less than 85 °C. Given a list of conditional terms (when, if, while) and a list of prescriptive terms (shall, must, should, have to, require, need, want, expect, wish or desire), we know the index of the conditional term and the prescriptive term by iterating through the tokens of a given requirement. Thus, to identify the pattern associated to a given requirement we use a set of rules that compares the index of the physical numeric dependencies, the index of the prescriptive term, and the index of the conditional term. This set of rules enables us to infer the Domain D and the Condition C.

A physical numerical value is sometimes followed by a tolerance. Thus, the patterns 1, 2 and 3 give rise to four more patterns where the domain D and the condition C are split into a nominal domain, a tolerance domain, a nominal condition and a tolerance conditione.g. pattern 5 follows from pattern 1.

Pattern 5: <Prescriptive> <Nominal Domain> <Tolerance Domain> <Nominal Condition> <Tolerance Condition> -PnDtDnCtC

The Control_Subsystem shall open the Inlet_Valve in 3 seconds +/-1 second, when the temperature of water in the Boiler is between 70 °C and 85 °C.

To make sure that two consecutive physical numerical dependencies form the so called <nominal, tolerance> pair of a domain D or a condition C, we check whether their units belong to the same physical dimension. For instance, in the requirement "When the temperature is less than 40°C, the pressure shall be less than Pa", the consecutives physical numerical dependencies <40 °C> and <30 Pa> do not belong to the same physical dimension since the former is a temperature (°C), whereas the latter is a pressure (Pa). However, in the requirement "The system shall control a pressure of 30 Mpa +/-5 Pa", the physical numerical dependencies <30 Mpa> and <5 Pa> belong to the same physical dimensiona pressure.

Finally, there are six more syntactic patterns according to whether there is a tolerance associated to the domain and/or the conditione.g. pattern 7 follows from pattern 1 and 5.

Pattern 7: <Prescriptive> <Domain> <Nominal Condition> <Tolerance Condition> -PDnCtC

The Control_Subsystem shall open the Inlet_Valve in less than 3 seconds, when the temperature of water in the Boiler is between 70 °C and 85 °C.

Step 12 < Tolerance Calculation>: The calculation of the minimum, maximum and nominal values defining the tolerance of a condition C or a domain D relies upon four patterns: (1) "X +/-Y" with X > Y, (2) "X more or less Y" with X > Y, (3) "from X to Y" with X < Y and (4) "between X and Y" with X < Y. If there is no tolerance, e.g. "the temperature shall be less than 50°C", the maximum and minimum values of the domain are identical. At present, there is a limit when the unit is not the same, e.g. "1 daN +/-10 N" or "from 10 Pa to 1 MPa", because we cannot compute the tolerance without using a unit convertor.

Step 13 <PBRs Modelling>: The NLP pipeline ends up with an XML specification that lists the PBRs in a structure that complies with the PPBRs data model in ENOVIA V6 Requirements. Thus, each PBR element has a statement, a nominal value, a minimal value and a maximum value that specify a domain D that can be inferred from the nominal domain and tolerance domain, a physical dimension and a unit attribute. The XML specification can finally be imported into ENOVIA V6 Requirements so as to automatically generate the PPBRs. This pure software development part of our proposal has not been implemented yet.

Once the PPBRs have been generated in ENOVIA V6 Requirements, designers can start the design synthesis thanks to the F, L and P layers of CATIA V6 RFLP.

Integrated, Parametric, Model-Based Product Design Synthesis

Design synthesis is the translation of input requirements into possible solutions satisfying those inputs [START_REF]INCOSE: Systems engineering handbook. A guide for system life cycle processes and activities[END_REF].

The design synthesis with the integrated CATIA V6 FLP product design environment consists in translating the input requirements (R) into functional, logical and physical solutions satisfying those inputs. In order to do so, designers recursively break down the functional requirements into functions (F) that transform flows. Functions are then implemented by dynamic logical units (L) that simulate the expected behaviour, whereas non-functional requirements that prescribe design constraints are implemented by inert structural organs (P). Once the design of a given hierarchical level is completed, we apportion the performance requirements of the current hierarchical level to the functions of the next lower level by either sticking with the same physical dimension (allocation) or by establishing new requirements resulting from specific implementation choices (derivation).

Parametric modelling enables designers to not only create flexible CAD model, but also to verify that the requirements comply with the design artefacts. When the PPBRs are directly specified in ENOVIA V6 Requirements, engineers have to manually create the requirements and the associated knowledge parameters. However, when requirements are imported from a document-based specification, the generation of PPBRs results from the NLP pipeline. The knowledge verification rules that link the PPBRs and the knowledge parameters of the design artefacts require domain-specific knowledge; consequently, they are manually defined by designers.

In the next section, we illustrate the transformation of TBRs into PBRs so as to generate PPBRs that drive the design synthesis of a mechanical ring. We use a mechanical ring as a case study because its design changes frequently and because it is a simple universally understood object.

3 Case Study 3.1 From Unstructured Specifications to PPBRs First, we put ourselves in the shoes of a contractor who wants to acquire a mechanical ring. We only write three requirements 10 to ease the illustration of our proposition. Each requirement is in a different specification (.doc, .pdf and SysML).

Then, we play the role of the OEM who receives the specifications. First, we upload the specifications and send them through the NLP processing pipeline. Once it has finished, the NLP processor generates the XML specification of PBRs.

In the XML specification, the data structure of the PBRs is defined in such a way that the PBRs are interpretable by ENOVIA V6 Requirements. Therefore, we can generate the PPBRs (Fig. 1) from the PBRs. Nevertheless, this functionality is not integrated into ENOVIA V6 because the NLP pipeline is a single capability of a broader on-going research project that requires Java EE as programming language.

Integrated, Parametric, Model-Based Product Design Synthesis

Now that the PPBRs are specified (Fig. 2(1)), designers can start the design synthesis of the ring with CATIA V6 RFLP. We only have non-functional PPBRs that prescribe design constraints, therefore they will be implemented by a structural design artefact, the ring, and more precisely, the external diameter, length and weight properties of the ring. By using parametric modelling, we define a design parameter for each property (Fig. 2(3)). A design parameter is a triple <Name, Physical dimension, Unit>e.g. a parameter named <D> whose physical dimension is <length> and unit is <mm> that drives the external diameter property of the ring.

After having associated the design artefacts' knowledge parameters with the geometrical features, designers define knowledge verification rules between the PPBRs and the design artefacts' knowledge parameters. For instance, Fig. 2(4) shows the knowledge verification rule verifying that the external diameter property of the ring belongs to the domain 20 mm +/-1 mm. It consists in defining a conjunction between two partial order relations "<" that constrain the design knowledge parameter <D> with the maximum and minimum values of the external diameter property (Fig. 2(4)) defined in the "Ring diameter" PPBR (Fig. 2(2)). As shows Fig. 2(5), a red light signals when the design does not comply with a requirement. In our case, the design artefact's knowledge parameter that stands for the external diameter property of the ring does not belong to the domain prescribed by the "Ring diameter" PPBR.

Results and Limitations

We conducted experiments with both real industrial and handcrafted data sets. The initial results that we obtained after analysing handcrafted specifications are encouraging. One key factor is that we wrote the requirements by following the requirements writing best practices defined in [START_REF]INCOSE: Guide for writing requirements. Requirements working group[END_REF]. Regarding the analysis of industrial specifications, the results are also promising, although we cannot fully validate the proposition without including a units converter.

Our solution presents a few limitations, such as the detection of sections when the headings functionality has not been used to edit .doc, .odf or .pdf. We were also challenged by original writing-style that came across while we were testing. The Stanford CoreNLP dependencies analyser relies upon statistics; consequently, it returns few false positive and false negative numerical dependencies. Finally, as previously explained, the translation of PBRs from TBRs is limited since the nominal and tolerance values must have the same unit.

Conclusion and Future Work

This paper presents a natural language processing pipeline to derive Parametric Property-Based Requirements from textual requirements.

In the future we plan to continue testing our algorithm on various specifications and integrate a unit converter to compute the minimum, maximum and nominal values defining the tolerance of a condition C or a domain D when the nominal and tolerance values do not have the same unite.g; 1 MPa +/-10 Pa. We will also develop the plug-in to load the XML specification into ENOVIA V6 Requirements so as to automatically generate the PPBRs from the PBRs.

  /CATIA V6 RFLP for Integrated, Model-Based Product Design

Pattern 2 :

 2 <Prescriptive> <Condition> <Domain> -PCD The Control_Subsystem shall, when the temperature of water in the Boiler is less than 85 °C, open the Inlet_Valve in less than 3 seconds. Pattern 3: <Condition> <Prescriptive> <Domain> -CPD When the temperature of water in the Boiler is less than 85 °C the Control_Subsystem shall open the Inlet_Valve in less than 3 seconds. Pattern 4: <Prescriptive> <Domain> -PD The Control_Subsystem shall open the Inlet_Valve in less than 3 seconds.

1 2Fig. 1 .

 11 Fig. 1. Two ENOVIA V6 Requirements windows: (1) the attributes to define a PPBR (left) and (2) the list of PPBRs that makes up the specification of the ring (right).

Fig. 2 .

 2 Fig. 2. (1) PPBRs, (2) external diameter PPBR, (3) design artefact's knowledge parameters, (4) knowledge verification rule, and (5) quantitative level of compliance.

https://www.modelica.org/.

http://www.3ds.com/products-services/catia/products/dymola.

http://www.3ds.com/fr/produits-et-services/catia/.

One should refer to[START_REF] Manning | Foundations of statistical natural language processing[END_REF] for further details on statistical natural language processing.

https://tika.apache.org/.

https://poi.apache.org/.

http://nlp.stanford.edu/software/tokenizer.shtml.

http://nlp.stanford.edu/software/tagger.shtml.

http://nlp.stanford.edu/software/stanford-dependencies.shtml.

(Req 1.) The diameter of the Ring shall be 20 mm +/-1 mm. (Req. 2) The length of the Ring shall be between 35 mm and 37 mm. (Req. 3) The Ring shall weight less than 500 g.