Gael Glorian
email: gael.glorian@labri.fr

Jean-Marie Lagniez
email: jean.marie.lagniez@huawei.com

Christophe Lecoutre
email: lecoutre@cril.fr

NACRE -A Nogood And Clause Reasoning Engine

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Many research efforts have been devoted to the design of general tools and algorithms for solving CSP (Constraint Satisfaction Problem) instances [START_REF]Handbook of Constraint Programming[END_REF]. The current mainstream approaches rely on complete methods that embed filtering procedures, search heuristics and learning mechanisms. Implementing a new technique into a CSP solver usually requires a strong knowledge of its whole structure. Moreover, to be fully accepted by the community, a comparison with the state-of-the-art is necessary, with careful re-implementation.

During the last two decades, the impressive progress in the related problem SAT (Satisfiability Testing), where variables are Boolean, has been achieved using nogood recording (called clause learning when nogoods are managed à la SAT) under a restart policy enhanced by a very efficient lazy data structure [START_REF] Moskewicz | Chaff: Engineering an efficient SAT solver[END_REF]. The value of clause learning has risen with the availability of large instances (encoding practical applications). Learning in SAT is an example of a successful technique derived from cross-fertilization between CSP and SAT: nogood recording [START_REF] Dechter | Enhancement schemes for constraint processing: Backjumping, learning, and cutset decomposition[END_REF][START_REF] Schiex | Nogood recording for static and dynamic constraint satisfaction problems[END_REF] and conflict-directed back jumping [START_REF] Prosser | Hybrid algorithms for the constraint satisfaction problem[END_REF] were originally introduced for CSP and were later imported into SAT solvers [START_REF] Bayardo | Using CSP look-back techniques to solve real-world SAT instances[END_REF][START_REF] João | GRASP: A search algorithm for propositional satisfiability[END_REF]. Progress in SAT has stimulated a renewed interest of the CSP community in nogood recording [START_REF] Boutaleb | no)good recording and robdds for solving structured (V)CSPs[END_REF][START_REF] Katsirelos | Generalized nogoods in CSPs[END_REF][START_REF] Lecoutre | Nogood recording from restarts[END_REF][START_REF] Richaud | Automata for nogood recording in constraint satisfaction problems[END_REF]. While some CSP solvers embed a clause learning component, some others require to use an external SAT solver or are built above a SAT solver. In some cases, the whole problem is encoded into CNF [START_REF] Soh | Proposal and evaluation of hybrid encoding of CSP to SAT integrating order and log encodings[END_REF][START_REF] Zhou | The Picat-SAT compiler[END_REF] or incrementally translated [START_REF] Feydy | Lazy clause generation reengineered[END_REF][START_REF] Ohrimenko | Propagation via lazy clause generation[END_REF]: the SAT solver often takes priority over the CSP solver for guiding search.

NACRE has been designed to be a hybrid solver, solving problem instances with either dedicated methods or SAT-inspired ones. Utilizing dedicated data structures, it can work efficiently with clauses (or nogoods); actually, the core of NACRE was conceived with this general objective in mind. Furthermore, NACRE has been conceived as a framework: it is designed to be upgradable, letting the user try new ideas on top of filtering or solving methods in a straightforward way. Moreover, it is possible to easily implement new reasoning methods concerning clauses and to specify minimization techniques to shorten them. Although NACRE has been conceived as a framework, its performances have not been left aside, as shown by the first place obtained in the CSP MiniTrack of the 2018 and 2019 XCSP3 competitions [START_REF] Boussemart | XCSP3: an integrated format for benchmarking combinatorial constrained problems[END_REF]. This track is for open-source mini-solvers1 and aims at discovering new ideas and solvers.

Preliminaries

A Constraint Network (CN) [START_REF] Lecoutre | Constraint networks: techniques and algorithms[END_REF] P is composed of a finite set of variables X and a finite set of constraints C. Each variable x has a (current) domain dom(x), which is the finite set of values that can be assigned to x; the initial domain of a variable x is denoted by dom init (x). Each constraint c involves an ordered set of variables scp(c), called the scope of c. A constraint c is semantically defined by a relation rel(c), which is the set of tuples allowed by (variables of) c. A tuple τ can be perceived as the instantiation of a subset of variables of X ; we note τ [x] the value of variable x in τ . A solution of a CN P is an instantiation of X such that all constraints c ∈ C are satisfied. A constraint c ∈ C is satisfied by an instantiation τ if and only if τ covers (i.e., instantiates) all the variables in the scope of c and τ corresponds to a tuple accepted by c. A CN is also usually referred to as a CSP instance.

To find solutions efficiently, it is important to filter the search space, typically by enforcing some local properties (called consistencies). Generalized Arc Consistency (GAC) [START_REF] Mackworth | Consistency in networks of relations[END_REF] is a very classical one. For a constraint c, it guarantees that each pair (x, a) with x ∈ scp(c) and a ∈ dom(x), admits a support on c, i.e., a tuple τ ∈ Π x∈scp(c) dom(x) such that τ [x] = a and τ ∈ rel(c). Constraint propagation is the process of enforcing a local consistency such as GAC, step by step, by soliciting all constraints to filter variable domains in turn until a fixed point is reached. A literal on a variable x corresponds to either an assignment of the form x ← a or a refutation of the form x a, with a ∈ dom init (x). The literals x ← a and x a are said to be positive and negative, respectively. In the following, we also denote literals with Greek letters (λ, δ, . . .). A decision δ is checked (in our algorithms) to be positive or negative by simply writing pos(δ) and neg(δ), respectively.

In this paper, we focus on backtrack search where decisions are taken in sequence to find solutions, and we assume that decisions always correspond to literals, which is a very classical approach in constraint programming. For example, MAC (Maintaining Arc Consistency) [START_REF] Sabin | Contradicting conventional wisdom in constraint satisfaction[END_REF] builds a search tree by systematically maintaining GAC after such (positive and negative) decisions. For the sake of precision, decisions refer to literals (unary constraints) added to the CN by the search algorithm when extending the current branch of the search tree, whereas deductions refer to inferences performed by the filtering process (e.g., the algorithm that enforces GAC). Deductions are also assumed to be literals (which is always the case with GAC).

Inside NACRE

We first succinctly introduce the data structures employed in NACRE before presenting the main components of the solver and the way nogood recording is conducted.

Data Structures

In NACRE, variables and values are represented through propositional variables: every pair (x, a) with x ∈ X and a ∈ dom(x) is paired to an integer representing the propositional variable x a (direct encoding [START_REF] Walsh | SAT v CSP[END_REF], where each CSP variable is represented by a sequence of propositional variables). Therefore, the variables are stored as virtual objects, meaning that each variable does not directly have its own domain but only knows an index where its domain starts (in a big contiguous table) as well as its initial domain size and its current domain bounds. There are two tables of this form (three in the case of clause learning). The first one, called varProps, stores the propositional variables and can be seen as the dense part of a full sparse set [START_REF] Aho | The Design and Analysis of Computer Algorithms[END_REF]. The propositional variables (referred to as PV) are objects that contain useful information such as the value of the PV into the original CSP variable, and also: (i) its current state (0 meaning that the value has been assigned to the CSP variable, 1 meaning that the value has been removed, and 2 when neither assigned nor removed), (ii) its local2 position in the table domValues described below, and (iii) a few Boolean flags to record if this is the last assigned value (by the solver), if it was used during the process of conflict analysis, etc. The second table, called domValues, is similar to a full reversible sparse set, i.e., it stores for each variable a two-part table separated by a limit which is the current size of the variable domain. The values on the left of this limit are still active into the variable domain, whereas the ones on the right have been removed. It allows us to remove and restore values easily and efficiently since it only requires two basic operations (swapping two integers and decrementing a limit) to either remove a value or assign a variable. This lightweight structure is almost backtrack-free as we just need to store the local position of the limit at each decision level.

Figure 1 shows an example on three variables x, y and z with dom(y) = {0, 1, 2} and dom(x) = dom(z) = {2, 4} as defined initially (initial state). The upper table is varProps, in which we can clearly see the values of the domain (we do not show the other PV information in the figure, for clarity). The lower table is domValues containing the local positions. We can observe that they are initially sorted in an increasing manner. Above these two tables, in light gray, we can see the integers associated to the propositional variables (which are essentially the indices of the table varProps).

Solver Components

Search Tree Exploration. NACRE provides various methods to solve CSP instances, given under XCSP3 format. Being a framework designed to handle clauses and nogoods, several stateof-the-art procedures are implemented. The first one is a standard MAC [START_REF] Sabin | Contradicting conventional wisdom in constraint satisfaction[END_REF][START_REF] Sabin | Contradicting conventional wisdom in constraint satisfaction[END_REF] approach, used with the method option -complete. It maintains arc-consistency while running a backtrack search based on binary branching. Binary branching means that at each search node, a left branch labeled with a positive decision x = a is developed first, and a right branch labeled with a negative decision x = a is developed next. To do this efficiently, NACRE uses time-stamps for both the constraints and the variables. Recall that a time-stamp is a value denoting the time at which a certain event occurred; time-stamps allow the progress of algorithms to be tracked over time. Technically, with time-stamps, a set is not required to represent the propagation queue, which renders propagation less expensive to deal with.

Algorithm 1: propagate(x: Variable): Boolean Q ← {x}; while Q = ∅ do y ← pickVariable(Q); forall c ∈ C | y ∈ scp(c) ∧ stamp[c] < stamp[y] do touched ← filter(c) ; if touched = WIPE OUT then return true; // Detected conflict Q ← Q ∪ touched; updateStamp(c); return false // No detected conflict
Algorithm 1 shows how propagation is conducted after the solver has decided on a variable x. Until a fixed point is reached, it picks (line 3) a variable y from Q (the propagation queue only contains x initially) by using Function pickVariable() that selects (and removes) a variable from Q. Then, we iterate over the constraints involving this variable (line 4), and the (time-)stamp of the current variable is compared with the (time-)stamp of the constraint. The rationale (discussed below) is that it allows us to determine whether it is useful to call the filtering algorithm (propagator), Function filter() at line 5, attached to the constraint c. This function returns the set of variables, touched, whose domains have been reduced by the propagator in order to add them to the propagation queue (line 7) or the special value WIPE OUT when a conflict occurs. Concerning stamping, note that the stamp of a constraint is updated at line 8, whereas the stamp of a variable is updated when it is added to Q (lines 1 and 7), meaning that its domain has just been modified (or touched ). This stamping system avoids unnecessary calls to constraint propagators, which could be very expensive. Indeed, when a constraint c is asked to filter, all variables in scp(c) are GAC regarding c (that is ensured by Function filter() at line 5) and stamp[c] is updated in such a way that stamp[c] ≥ stamp[x] for all x ∈ scp(c).

It is clear that it is unnecessary to consider c again before any variable of its scope is touched. When a variable x is touched, then stamp[x] is updated, and x is added to the propagation queue. Therefore, if a variable x of scp(c) is touched after c is revised then stamp[x] > stamp[c] and consequently the condition (line 4) is satisfied when x is picked from the queue and the constraint is considered again to be revised.

Search Heuristics. In NACRE, we have a small pool of variable ordering and value ordering heuristics for tuning search. The variable ordering heuristic dom [START_REF] Bitner | Backtrack programming techniques[END_REF][START_REF] Haralick | Increasing tree search efficiency for constraint satisfaction problems[END_REF] (option -dom) was the first popular dynamic heuristic. It uses the current size of the variable domains to select the next variable to instantiate; the smaller, the more preferred. A second variable ordering heuristic is dom/deg [START_REF] Bessière | MAC and combined heuristics: Two reasons to forsake FC (and CBJ?) on hard problems[END_REF] (option -domdeg), a dynamic heuristic that combines both the domain size and the degree of the variables. The degree of a variable is the number of constraints where the variable is involved. This heuristic selects the variable with the smallest ratio of its current domain size to its degree. A third variable ordering heuristic is the robust adaptive heuristic dom/wdeg [START_REF] Boussemart | Boosting systematic search by weighting constraints[END_REF][START_REF] Wattez | Refining constraint weighting[END_REF] (option -domwdeg, used by default) that aims at choosing in priority the variables involved in the most falsified constraints. Technically, in each constraint object, a counter (initialized to 1) is incremented whenever the constraint is at fault (i.e., fails while filtering) during the search.

Once a variable is selected, a value must be chosen from its domain. Different value ordering heuristics are available in NACRE. In addition to the classical value heuristics (min-value, maxvalue, . . . ), a weak version of the SAT-based heuristic polarity has been implemented [START_REF] Pipatsrisawat | A lightweight component caching scheme for satisfiability solvers[END_REF]. This heuristic, called saving option (-saving), tries to assign the selected variable to the last value assigned to it. If the choice is not valid anymore, meaning that the value has already been deleted, then the underlying value heuristic is solicited. Two main options are available, each one having its opposite. On the one hand, with min-value (option -valMin), we choose the minimal value in the domain of the variable; it is also possible to select the max-value inside NACRE (option -valMax). Choosing the minimal or maximal value, depending on the problem, can greatly improve the performance of the solver. On the other hand, the heuristic first-value (option -valFirst), called lexico sometimes, is a second way of selecting values, based on the implementation of the domains. As we use structures similar to sparse sets, we can select the first PV inside the domValues set. We can also select the last one before the limit (option -valLast). It is also possible to fully randomize value selection (option -valRand).

Restart Policies. NACRE has some built-in restart policies that allow us to test nogood or clause learning under different settings. Because restart policies can have a strong impact on search efficiency, we have implemented the most popular policies, from geometric progression to Luby sequence [START_REF] Luby | Optimal speedup of las vegas algorithms[END_REF]. The first restart policy available is the no restarting policy (option -noRst), sometimes useful to be used as a baseline. Then comes the Luby-based sequences: first terms are of the form (1, 1, 2, 1, 1, 2, 4, 1, . . .). There are 3 available Luby-based sequences, Luby(#rst) × N where N can take the values 10, 50 and 100 (options are respectively -luby10, -luby50 and -luby100) and #rst is the number of runs already done in order to get the correct term from the sequence. This is useful for trying more or less aggressive restarts policies. Finally, geometric sequences are available for a more constant progression of the cutoff value. Four of them are implemented inside NACRE. The default one has a 3% progression, others are 10% (option -10perc), 50% (option -50perc), 100% (option -double). It allows the user to choose from a wide variety of policies, even possibly changing it during search with a few lines of code.

Nogood Recording

The benefit of recording nogoods [START_REF] Schiex | Nogood recording for static and dynamic constraint satisfaction problems[END_REF] is to avoid some form of thrashing [START_REF] Gomes | Heavy-tailed phenomena in satisfiability and constraint satisfaction problems[END_REF], i.e., exploring the same unsatisfiable subtrees several times. There are two classical methods to identify and store nogoods: during the search or at restarts. Nogood from restart extracts information by analyzing the current search tree when a restart is triggered. Two well-known and effective methods have been implemented into our framework, the negative last-decision nogoods [START_REF] Lecoutre | Recording and minimizing nogoods from restarts[END_REF] (nld-nogoods) and the increasing-nogoods [START_REF] Jimmy | Increasing nogoods in restart-based search[END_REF]. Nld-nogoods (option -nld) are extracted at each restarts from the last branch of the search tree. Let us assume that the sequence of labels all along the rightmost branch of the current search tree being developed is Σ = δ 1 , . . . , δ m , where each decision of Σ is either a positive or a negative decision. It is known that for any i such that 1 ≤ i ≤ m and neg(δ i ), the set ∆ = {δ j : 1 ≤ j < i ∧ pos(δ j )} ∪ {¬δ i } is a reduced nld-nogood. A reduced nld-nogood {δ i , δ j } can be written δi∈∆ δ i ⇒ ¬δ j in directed form. Increasing-nogoods (method option -incng) are a compressed form of all the reduced nld-nogoods that can be extracted at each restart, without losing any effectiveness. Indeed, it is rather easy to deduce that there are some similarities between nld-nogoods: the left part of the implication is shared between several nogoods. Then, it is possible to rewrite the set of nld-nogoods as a constraint and use a dedicated propagator that ensures GAC.

Another way to learn nogoods is during search when a conflict occurs. The concepts of nogoods and clauses are closely related. The former is mostly used in the CSP context and the latter in the SAT one. NACRE implements a generic clause reasoning engine using lazy explanations (option -ca). We consider the g-nogood learning procedure [START_REF] Katsirelos | Unrestricted nogood recording in CSP search[END_REF], and the lazy explanations, proposed in [START_REF] Gent | Learning when to use lazy learning in constraint solving[END_REF][START_REF] Gent | Lazy explanations for constraint propagators[END_REF] as a starting point. A way of generating a nontrivial nogood at each conflict consists in deeply analyzing the sequence of propagation steps by constructing an implication graph à la SAT and identifying a reason of the failure [START_REF] Zhang | Efficient conflict driven learning in boolean satisfiability solver[END_REF]. An implication graph is a directed acyclic graph (DAG) that records the relationships existing between literals, as can be observed during the solving process. Each vertex represents a literal λ (with its associated level), and its incoming arcs in the DAG represent the reasons that force it. When a conflict occurs, each cut in the implication graph that leads to the conflict can explain it. This set of literals, used as a conflict explanation, can be blocked in the future by deriving a clause from it. Generally, we are only interested in Unique Implication Points (UIPs) [START_REF] João | GRASP: A search algorithm for propositional satisfiability[END_REF] that are vertices at the current decision level dominating the conflict. Importantly, they can be used to safely perform non-chronological backtracking (a form of back jumps) up to the decision level that is the maximum value among all decision levels associated with the literals of the conflict explanation.

In practice, the explicit construction of the implication graph can be avoided. Indeed, conflict clauses can be derived from the current formula if we simply save independently the explanation expl(λ) of each deduction λ. More precisely, assuming that all such explanations are stored, and starting from a conflicting set Λ of literals, the conflict analysis procedure iteratively selects a literal λ ∈ Λ implied at the most recent decision level (i.e., the current largest decision level) and replaces λ by expl(λ) in Λ. This process is stopped when only one literal from the current level remains in Λ. Consequently, the literals remaining in Λ, after the procedure is complete, forms a UIP. Because Λ represents an explanation for the conflict, the clause λ∈Λ ¬λ is a logical consequence of the problem (at least one literal of Λ must be false). We can observe that upon backtracking this clause remains unit (with literal λ unit ) until the highest decision level of the other literals in Λ \ {λ unit } is reached.

On one hand, SAT formalism is very simple (each time a literal is deduced, its explanation can be computed by considering the clause that triggers it) and very well suited to nogood extraction by an analysis of the implication graph. On the other hand, CSP formalism is more sophisticated (e.g., dealing with global constraints) and makes more difficult the precise identification of the reason of a simple deduction: building the implication graph is no more straightforward. Moreover, even if we may consider constructing such graphs, computing the explanation of each deduction can drastically slow down the search algorithm if no proper care is taken. In NACRE, a generic non-intrusive method that allows us to compute explanations has been implemented.

Explanations are generated on demand in three situations. Firstly, if a literal λ is propagated from a clause Λ, then the explanation is given by the negation of literals of Λ \ {λ}. Secondly, a literal can be propagated by the domain constraints: atMostOne and atLeastOne. The atLeastOne constraint is used as explanation when only one value a is remaining is the domain of x, in this case we have expl(x = a) = {x = b : b ∈ dom init (x) \ {a}}. The atMostOne constraint is used as explanation when a value b is deleted from x when x is assigned to a, in this case we have expl(x = b) = {x = a}. Finally, a literal λ from x can be propagated while running the filtering algorithm φ c associated with a constraint c that is not a clause. In this case, it is possible to explain why λ has been propagated by considering the deleted values for the variables of scp(c) \ {x}. To avoid systematic computation of explanations at each propagation, NACRE uses a stamp system that associates an integer with each value. When a value a from x is deleted its stamp is updated in such a way that stamp(x, a) is greater than all stamp(x , a ) such that a is deleted from x . In order to be backtrack free, we use an integer that is incremented at each deletion. Thus, when an explanation is required for a value a from x, it is enough to keep in memory which constraint has been used to perform the propagation. Indeed, it is enough to consider the value a of x ∈ scp(c) \ {x} such that stamp(x, a) > stamp(x , a ). Once the conflict analysis process is over, a new clause is added to the nogood base and a non-chronological backtrack is performed. Of course, for managing the nogood base, we adopt the classical two-watched-literals strategy [START_REF] Moskewicz | Chaff: Engineering an efficient SAT solver[END_REF].

Algorithm 2 is an variation of Algorithm 1 that handles clause propagation too. We choose to fully propagate the constraints (line 3 to 10) before the clauses (line 12 to 18). Indeed, playing with clauses and nogoods can lure an adaptive heuristic like dom/wdeg [START_REF] Glorian | Combining nogoods in restart-based search[END_REF]. Triggering in priority conflicts from the CSP core avoids this kind of behavior. These two propagation processes are done until a fixed point is reached (when clause propagation loopbcp-becomes useless). The constraint propagation loop is similar to Algorithm 1 with one slight difference at line 5. The picked variable is added into a new queue Q SAT . It allows the loop at line 12 to iterate over all the variables that were touched during GAC enforcement.

Clause Database Reduction. The procedure that ensures unit propagation on nogoods (clauses) highly depends on the size of the base. To maintain manageable the base of nogoods, and perform unit propagation at a reasonable cost, it is common to reduce the base by deleting clauses considered to be irrelevant to the next search steps. Several measures have been proposed for this purpose [START_REF] Audemard | On freezing and reactivating learnt clauses[END_REF][START_REF] Audemard | Predicting learnt clauses quality in modern SAT solvers[END_REF][START_REF] Eén | An extensible sat-solver[END_REF]. Generally, SAT solvers use one of the following strategies: activity that considers a learned clause as irrelevant if its activity or its involvement in recent conflict analysis is marginal [START_REF] Eén | An extensible sat-solver[END_REF]; or literal block distance (LBD) that uses the number of different levels involved in a given learned clause to quantify the quality of learned clauses [START_REF] Audemard | Predicting learnt clauses quality in modern SAT solvers[END_REF]; clauses with smaller LBD are considered more relevant.

Even if these measures are well adapted for SAT, there are not well suited in our case: a hybrid solving method requires a hybrid clause quality measure. Indeed, it is often the case that clauses contain a large number of literals. But contrary to SAT, it is possible to identify literals that are linked by a common CSP variable. In this case, clauses that are related to few CSP variables are more powerful regarding unit propagation. To take this information into account, we propose to use a new measure that considers the number of CSP variables involved in the clause through its associated propositional variables. When two clauses have the same number of CSP variables, then their activity is used to decide between them. An important point regarding the reduction of the base concerns its frequency. In NACRE, we chose to set the initial size limit of the base to 4, 000 nogoods (pragmatic choice, made from experiments). When the limit is reached, half of the learned clauses are removed using our custom measure.

Figure 1 :

 1 Figure 1: Data Structures for Variables/Values in NACRE.

NACRE source is available at https://github.com/crillab/nacre_mini.

We call local a position that has been normalized according to the current variable domain starting index.

The -cm option stands for verbose mode 0: only the result and the solution, if one is found, are displayed.

PicatSAT was not an official competitor but is given as a comparison with the main competition track.

Generated from the raw data of the official competition website: http://www.xcsp.org/competition.

This work has been supported by the project CPER Data from the region "Hauts-de-France".

Algorithm 2: propagate hybrid (x: Variable): Constraint while

until ¬continue; return null ; // no detected conflict

After each reduction, the size limit is increased by 500.

NACRE at XCSP3 Competitions

NACRE was submitted to the CSP MiniTrack of the XCSP3 competition [START_REF] Lecoutre | Proceedings of the 2018 XCSP3 competition[END_REF]. This track is designed for sequential open-source software and allows constraint solver developers to enter the competition without having to implement all constraints from the so-called XCSP3-core. Solvers are evaluated on a restricted set of constraints: intension, extension, allDifferent, sum and element. Solvers were run on a cluster of Xeon@2.67GHz with 32GiB of memory. Sequential solvers were allocated 15500 MiB of memory and a time limit of 2400 seconds. The 2018 (resp. 2019) CSP MiniTrack was a selection of 176 (resp. 200) instances. NACRE competed in its hybrid form with the following options: -ca -luby100 and -cm 3 . It enables our conflict analysis method and uses the Luby sequence as restart policy. In 2019, we also submitted the standard MAC version (options: -complete -l100 -cm) to the competition in order to compare them.

Table 1 5 shows that, in 2018 (resp. 2019), NACRE solved 49% (68%) of the 176 (200) instances, which represents 76% (78%) of the Virtual Best Solver (VBS). It means that 76% (78%) of the instances that can be solved by any competing solvers are solved by NACRE. We can also see that NACRE performs well on both SAT and UNSAT instances. The impact of clause learning is stronger on the 2019 results (17% more solved instances).

Conclusion

In this paper, we have presented NACRE, a nogood and clause reasoning engine. Our main goal is to provide the community with an expansible tool, useful in particular for easily experimenting with nogoods and clauses. We paid attention to practical performances, and we demonstrated that in the 2018 and 2019 XCSP3 competitions. For the next major version, we project to handle all constraints from XCSP3-core.


Dr Arthur Mageau
email: arthur.mageau@inserm.fr

MD Geoffrey Aldebert

Damien Van Gysel

MD Thomas Papo

MD Jean-François Timsit

PhD Karim Sacré

Sars-Cov2 infection among Systemic Lupus Erythematosus inpatients in France: a nationwide epidemiologic study

Since the global emergence of Sars-Cov2 at the end of 2019, a special concern has raised regarding patients with rheumatic and inflammatory diseases, such as systemic lupus erythematosus (SLE) [START_REF] Contributors | /CRI/IMIDIATE consortium and. Severity of COVID-19 and survival in patients with rheumatic and inflammatory diseases: data from the French RMD COVID-19 cohort of 694 patients[END_REF]. Indeed, many treated SLE patients are immunocompromised and often suffer from chronic kidney or cardiovascular diseases [START_REF] Mageau | Septic shock among patients with systemic lupus erythematosus: Short and long-term outcome. Analysis of a French nationwide database[END_REF]. Recently, anti-interferon antibodies have been implicated in severe Sars-Cov2 infection [START_REF] Bastard | Autoantibodies against type I IFNs in patients with life-threatening COVID-19[END_REF] while it has been known for decades that SLE patients may produce such autoantibodies [START_REF] Gupta | Association between anti-interferon-alpha autoantibodies and COVID-19 in systemic lupus erythematosus[END_REF]. Although available data from short SLE series are reassuring [START_REF] Mathian | Clinical course of coronavirus disease 2019 (COVID-19) in a series of 17 patients with systemic lupus erythematosus under long-term treatment with hydroxychloroquine[END_REF][START_REF] Ramirez | COVID-19 in systemic lupus erythematosus: Data from a survey on 417 patients[END_REF], large-scale nation-wide studies are still needed to assess the risk of developing severe Sars-Cov2 infection in SLE.

We therefore used the French healthcare database system called "Programme de Médicalisation des Systèmes d'Information" (PMSI) -which contains hospitalization data of all inpatients in France -to analyse the SLE population that had at least one stay in a French hospital between March and October 30 th , 2020 (Supplementary Material). On this population, we compared SLE patients with or without Sars-Cov2 infection (SLE/COVID-19 + , SLE/COVID-19 -). Among SLE/COVID-19 + , we distinguished patients with poor outcome after Sars-Cov2 infection and patients with good outcome after COVID-19. We defined poor outcome as admission to Intensive Care Unit (ICU) or death. We also compared the in-hospital mortality associated with Sars-Cov2 infection in SLE and in the total population in France.

Based on the 10th International Classification of Diseases "M32" and "L93" diagnosis code, we identified 11 055 SLE patients who had at least one stay in a French hospital between March 1 st and October 31 st , 2020 (Supplementary materials, Figure S1). Among them, 1 411 (12.8%) also had a COVID-19 diagnosis code. Characteristics of SLE/COVID-19 + and SLE/COVID-19 -patients are given Table 1. These 1 411 SLE/COVID-19 + patients experienced 1 721 hospital stays during the period of study. Among these SLE/COVID-19 + stays, 293 (17 %) took place in Intensive Care Unit (ICU).

The mean Simplified Acute Physiology Score (SAPS) II at admission was 35.4+/-16.8. In ICU, 78 (26.7%) and 71 (24.7%) SLE/COVID-19 + patients required invasive or non-invasive mechanical ventilation, respectively. Overall, 134 (9.5 %) SLE patients admitted for COVID-19 died. The in-hospital mortality rate was almost 4 times higher in SLE/COVID-19 + as compared to SLE/COVID-19 -patients admitted during the same period (9.5 % vs 2.4 %, p <0.001). Interestingly, while the overall mortality rate was lower in SLE/COVID-19 + patients as compared to the total population admitted for Sars-Cov2 infection in France during the same period (9.5 % vs 15.7%, p <0.0001), the mortality rate at a younger age appeared higher in SLE patients after stratification according to the age (Supplementary materials, Figure S2).

Our study based on a comprehensive nationwide database confirms that patients with SLE are more likely to develop severe Sars-Cov2 infection when they have comorbidities already identified as risk factors of severe infection in the general population, such as older age, male gender, and hypertension [START_REF] Contributors | /CRI/IMIDIATE consortium and. Severity of COVID-19 and survival in patients with rheumatic and inflammatory diseases: data from the French RMD COVID-19 cohort of 694 patients[END_REF]. Poor outcome was associated with chronic kidney disease (CKD) status thus confirming that CKD increases the risk for severe infection in SLE and is a major predictor of mortality and morbidity in these patients [START_REF] Mageau | Septic shock among patients with systemic lupus erythematosus: Short and long-term outcome. Analysis of a French nationwide database[END_REF].

Given the importance of male sex as a poor prognosis factor of COVID-19, the lower mortality rate observed among SLE population may be explained by an unbalanced sex ratio (F/M 8.5:1.5). On the other hand, the mortality of SLE/COVID-19 + patients seems higher in the youngest patients as compared to the general population with Sars-Cov2 infection.

Because lupus activity and the need for immunosuppressive drugs decline with age, the higher mortality rate observed in younger patients suggests that SLE disease impact COVID-19 outcome. This finding argues for prioritizing vaccination studies and policies in the youngest patients with SLE.

AM designed and conducted analysis and wrote the manuscript. GA, NC, DG and TP were involved in the project development and edited the manuscript. KS and JFT directed the project and wrote the manuscript.

PATIENT AND PUBLIC INVOLVEMENT

Patients or the public were not involved in the design, or conduct, or reporting, or dissemination plans of our research

TABLE Table 1 Characteristics of SLE inpatients in France with and without COVID-19 infection between March and October 2020.

 Table Level of significance for the difference between SLE/COVID-19 + with good outcome vs SLE/COVID-19 + with poor outcome. Student's T test was performed for quantitative variables and chi-square test was used for categorical variables.

		SLE/COVID-19 -	SLE/COVID-19 +
			good outcome	poor outcome	p*
		n= 9644	n= 1070	n= 341
	Age, years (mean +/-sd)	52.1+/-17.2	55.8+/-18.7 61.5+/-17.6	<0.0001
	Male gender, n (%)	1 144 (16 .1)	171 (16.0)	96 (28.2)	<0.0001
	Chronic kidney disease, n (%)	1312 (13.6)	184 (17.2)	108 (31.7)	<0.0001
	Diabetes, n (%)	1251 (13.0)	173 (16.2)	72 (21.1)	0.03
	High blood pressure, n (%)	3404 (35.3)	480 (44.9)	215 (63.1)	<0.0001
	Obesity, n (%)	1811(18.8)	240 (22.4)	86 (25.2)	0.29
	Chronic pulmonary disease, n (%)	1275 (13.2)	240 (22.4)	102 (29.9)	0.005
	Cardiovascular event history, n (%)	2050 (22.8)	291 (28.5)	140 (43.8)	<0.0001
	Lupus nephritis history, n (%)	1271 (13.2)	158 (14.8)	77 (22.6)	0.0007
	Sjögren syndrome, n (%)	1088 (11.3)	124 (11.6)	38 (11.1)	0.8
	Antiphospholipid syndrome, n (%)	803 (8.3)	107 (10.0)	144 (12.9)	0.13
	Solid organ transplantation, n(%) **	212 (2.2)	22 (2.1)	14 (4.1)	0.04

* ** Solid organ transplantation during the last 10 years.

FUNDING

PhD fellowship support for Arthur Mageau was provided by Agence Nationale pour la recherche (n°ANR-19-CE17-0029)

AUTHOR CONTRIBUTIONS


Alain Sellier
email: alain.sellier@insa-toulouse.fr

Thierry Vidal
email: bvidal@insa-toulouse.fr

Hugo Gérard

Jacques Cagnon

Laurie Buffo-Lacarriere
email: buffo-lacarriere@insa-toulouse.fr

Stéphane Multon
email: estephane.multon@insa-toulouse.fr

Hugo Cagnon
email: chugo.cagnon@insa-toulouse.fr

Requirements for the Modeling of Medium-Term Behavior of Nuclear Containment Concrete for a "Loss of Coolant Accident" Analysis

Keywords: nuclear containment, concrete, creep, temperature, AAR, DEF, cracking

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

Introduction

The French Nuclear industry provides 79% of the electricity used in France, but it has to face with two mains problems: first the aging of its nuclear power park which, for some nuclear power plants, arrive at their design lifetime; and secondly there is an increasing need of safety demonstration following the severe accidents of Chernobyl and Fukushima. So, the different nuclear agencies, and the French government work together tirelessly to control the nuclear risk. Among the new measures already begun, there are several research programs aiming at increasing the knowledge level to predict more and more accurately the behavior of a double wall concrete vessel in case of loss of coolant accident (LOCA). If such an accident would occur, the vapor would be stopped by the double concrete wall system, and redirected to filters which work under a 4 bars relative pressure. At this pressure, the vapor temperature inside the vessel can reach 180°C, and be maintained until two weeks. This vapor pressure and this temperature constitute the internal conditions in the scenario of LOCA. The current research programs MACENA and ENDE, funded by the French research agency (ANR) work on tools to predict accurately the double wall mechanical behavior in case of LOCA. The objective of the research project MACENA is to provide an accurate finite element model able to predict if a vapor leakage would be possible in such conditions. The fact that some nuclear power plants have more than 50 years old has to be considered in the models, particularly because the internal concrete wall is pre-stressed with steel wires. With time, their pre-stress are relaxed, and on the other hand, the concrete, constantly compressed by the wires, creeps in compression, leading to a concrete shrinkage which contributes also to the loss of pre-stressing. While the pre-stressing is able to balance the tensile stress induced by the internal pressure in case of LOCA, the leak risk is low, but, if due to the loss of pre-stressing, the concrete passes in tension, even locally, there is a risk of cracking, and then a risk of radioactive vapor leakage. Of course, the concrete can itself resist to a tensile stress of a few MPa, depending on its own tensile strength; this strength constitutes then a supplement safety margin which can be considered for the leakage prediction in case of LOCA. But the current tensile strength of a given concrete 50 years old is not so easy to know for different reasons: first, the tensile strength of concrete is strain dependent; it decreases if the strain exceeds the tensile strain, leading to a progressive strain localization causing finally a localized crack. So, if earlier in the wall life, a large tensile strain has been reached, even during a very short period (for example at early age before the pre-stressing, or during a decennial pressurization test), the tensile strength is definitively reduced. Another problem resides in the heterogeneity of concrete, which leads to an intrinsic variability of the tensile strength impossible to neglect in case of large structures as shown during the CEOS.fr research program [START_REF]French National research project for design and assessment of special concrete structure toward cracking and shrinkage[END_REF]. A last problem, already observed in North America, concerns the concrete chemistry: in fact the cement paste and the aggregates of concrete can evolve chemically since the vessel construction. Among the different possibilities, the alkali aggregate reaction, the delayed ettringite formation, the corrosion of steel bars due to ingress of carbonate or other chemical species, could occur. Finally, the prediction of the double wall containment vessel behavior under LOCA must be based on a numerical simulation considering the whole causes of pre-stressing loss and tensile strength reduction previous to the envisioned LOCA. The paper is a listing attempt of model characteristics needed in order to predict the cracking under LOCA conditions.

General principle of a model

As mentioned in the introduction, several phenomena influencing the cracking occurrence have to be considered to predict the risk of radioactive vapor leak in case of LOCA. They can be divided in three main sets: Loading effects prior to the LOCA, loading effects during the LOCA, and material characteristics, as illustrated in Fig. 1. Even if each of these phenomena can occur independently, their consequences in terms of mechanical damage, inelastic strains, chemical evolution are most of the time irreversible. So, the material properties just before the LOCA depend on the solicitation history since the concrete pouring. Several computational strategies can be adopted to consider the damage accumulation; for example several models can be chained, a first one dedicated to the prediction of cracking at early age, a second one taking into account concrete creep and wires relaxation, and a last one able to consider the effects of temperature and vapor pressure to assess the cracking occurrence during the LOCA. Therefore, chaining several models leads rapidly to complications in terms of continuities of material state between two consecutive modeling, particularly if the internal state variables are not exactly the same between two consecutives models (for example the damage variable can be a tensor in a model and a scalar in another one, the permanent strains can be split in a different number of tensors depending on the model used). To avoid these problems of links between elementary models, a better strategy consists to consider all the afore-mentioned phenomena relative to a given material in a single constitutive model. The main advantage is the continuity, along the numerical simulation, of the internal state variables, which considers then automatically the damage accumulation until the LOCA. The following sections focus on the constitutive model clarified in LMDC Toulouse (France) to consider the different phenomena involved in the concrete behavior since its pouring until the LOCA. Three main phenomena will be discussed: (1) the creep, (2) the way to consider the influence of some chemical reactions such as hydration or alkali aggregate reaction on the mechanical behavior law, (3) the probabilistic scale effect. The paper ends with considerations about coupling between these phenomena.

The concrete creep modeling requirements

Multiaxiality and logarithm form of the long term creep. The internal wall of the structure is pre-stressed in two directions, with a compression stress close to 12 MPa and 8MPa in horizontal and vertical direction; so the creep model has to be able to manage the multi-axiality of the creep. Another particularity is the long duration of loading. Some nuclear power plants are around 50 years old, and their creep deformations are not yet stabilized, presenting all the same a slow decrease of the creep velocity with the time. The FIB recommendations, following the code model 2010 [START_REF] Walvaren | Model Code[END_REF], proposes using a logarithm function of the time to model this type of non-asymptotic creep; it is convenient to compute concrete elements under quasi-constant loading, but not really adapted for a numerical implementation for which a differential formulation with internal variables is more suitable. Indeed, a differential formulation is able to consider transient changes of loading, temperature or humidity without any problem, while a creep function would be added for each change weighting the formulation with the complexity of the loading.

( M ij c M E ij M ij C t ε τ ε ε = ∂ ∂ (1)
In this model, the internal variables are the directional consolidation coefficients ( )

( M ij c C ε in Eq.1)
which evolve with the creep strain main values, affecting the apparent viscosity (i.e. the characteristic time in Eq.1) of a Maxwell module (Fig. 2). The initial anisotropy of the creep is induced by the use of the elastic strain E ij ε in Eq.1, while the induced anisotropy is the consequence of the coefficients of consolidation ( )

( M ij c
C ε ) evolution (see ref [START_REF] Sellier | Concrete creep modelling for structural applications : non-linearity, multi-axiality, hydration, temperature and drying effects[END_REF] for more details).

Effect of humidity and temperature on creep. Concerning the influence of temperature and humidity on creep velocity, it is necessary to consider them properly. In fact, as illustrated in [START_REF] Cagnon | Drying creep in cyclic humidity conditions[END_REF] and reported in Fig. 3, the relative humidity variations induces drying creep during drying periods, and swelling during wet periods; more, the creep velocity depends on the moisture content. The LMDC creep model considers the effect of humidity through the saturation degree which reduces the creep velocity, while the drying creep and the drying shrinkage are considered thanks to poro-mechanics considerations. Indeed, with a poro-mechanic formulation the effects of drying or wetting are the consequences of the interstitial pressure variations (Fig. 2).

Effects of chemical reactions on concrete mechanical behavior

The chemical reactions can be divided in two main types, the ones affecting chemically the solid skeleton of concrete such as hydration or leaching [START_REF] Buffo-Lacarrière | Finite element modelling of hardening concrete: application to the prediction of early age cracking for massive reinforced structures[END_REF], and the ones constituting an internal loading through the pressure induced into the pores by the formation of new components such as gel of alkali aggregate reaction [START_REF] Multon | Cement and Concrete Research Multi-scale analysis of alkalisilica reaction ( ASR ): Impact of alkali leaching on scale effects affecting expansion tests Cem[END_REF] or delayed ettringite [START_REF] Salgues | DEF modelling based on thermodynamic equilibria and ionic transfers for structural analysis[END_REF]. Chemical reactions affecting the solid skeleton: as explained in [START_REF] Buffo-Lacarrière | Finite element modelling of hardening concrete: application to the prediction of early age cracking for massive reinforced structures[END_REF], both the precipitation or the dissolution of hydrates can be considered with a differential form of the behavior law, but the term corresponding to the stiffness variation ( t S ijkl ∂ ∂ / in Eq.2) must be considered only for leaching or dehydration (i.e. if the Young modulus (E) variation is negative in Eq.2). It corresponds to the loss of elastic energy associated to the dissolution of a loaded hydrate.

E kl ijkl E kl ijkl ij t S t E H t S t ε ε σ ∂ ∂       ∂ ∂ - + ∂ ∂ = ∂ ∂ ' ~(2)
In Eq.2, ' ~ij σ is the effective stress (the stress in the solid skeleton of concrete), E kl ε the elastic strain and H the Heaviside function. The evolution of mechanical parameters such as Young's modulus, Poisson ratio or strengths can be considered thanks to a relationship between the hydration degree and the corresponding parameters. A useful form for this type of relationship is the De-Shutter one (Eq.3) [START_REF] Schutter | Degree of hydration-based description of mechanical properties of early age concrete[END_REF].

n ref ref ref E E         - - = ξ ξ ξ ξ 1 ) ( (3)
Eq. 3 is presented here for the Young's modulus, but similar relations exist for all the mechanical parameters. In these relations, ref ξ is the solidification threshold, constant whatever the parameter for a given concrete, and n the non-linearity exponent, different for each parameter but quasiindependent of the concrete formulation. Chemical reaction with pore pressures: Alkali reaction and delayed ettringite formation belong to this second family of chemical reaction affecting the concrete behavior. The coupling between these reactions and the other aspects of the concrete behavior can be advantageously formulated in the framework of the poro-mechanics, in order to benefit from other features of the model; like creep or hydration-dehydration effects which are then "encapsulated" in the solid skeleton behavior law, while the AAR or DEF are added as internal swelling pressure acting in the concrete pores. A typical form for the poro-mechanics formulation is given by Eq.4.

∑ = ∂ ∂ - ∂ ∂ = ∂ ∂ N i i i ij ij t P S B t t 1 ) ( ' σ σ (4)
with B the Biot coefficient of the porous media, N the number of phases in the porosity, S i the saturation degree of the porosity with the phase i, P i the internal pressure induced by a neo-formed phase (AAR gel or DEF), the water and the gas pressure. In [START_REF] Grimal | Creep, Shrinkage, and Anisotropic Damage in Alkali-Aggregate Reaction Swelling Mechanism-Part I: A Constitutive Model ACI Mater[END_REF][START_REF] Grimal | Concrete modelling for expertise of structures affected by alkali aggregate reaction[END_REF], the form adopted to assess the pressure induced by AAR gel and DEF is Eq.5:

( ) ( )                 ∂ ∂ + ∂ - ∂ + ∂ ∂ - ∂ ∂ = ∂ ∂ t tr S t tr S B t P P t A M t P i pl i i pl i i i i i i i i , , lim 0 . . . ε ε ε φ φ (5)
In Eq.5 , M i is the Biot modulus for the phase i, i φ the final amount of the chemical reaction product, i 0 φ is the volume not efficient for the pressure (gel in pore connected to reactive sites, or gel with insufficient molar volume [START_REF] Bouzabata | Swellings due to alkali-silica reaction and delayed ettringite formation: Characterisation of expansion isotropy and effect of moisture conditions[END_REF]), A i is the reaction advancement ranging from 0 at the beginning of the reaction until 1 when the reaction is ended. i P lim is the limit pressure before cracking in free swelling conditions. i pl, ε is the plastic strain corresponding to the reaction products solidified in cracks caused by an excessive pressure. This plastic strain is computed thanks to a plastic criterion written in terms of effective stress concentration in the vicinity of the reactive sites:

i i I i i I i g Rt Rt P k f - + = ) , min( , σ (6)
In [START_REF] Vidal | Effect of Temperature on the Basic Creep of High-Performance Concretes Heated between 20 and 80[END_REF], k i is the stress concentration coefficient, ) ( ,i pl i Rt ε the effective tensile strength of the solid skeleton. The function 'min' is used to consider only admissible effective stresses in the skeleton since another criterion is used to limit this last. Coupling between creep and chemical reactions: finally, the coupling between chemical reactions and other phenomena is automatically considered since the effective stress used in Eq.6 takes into account Eq. 2 and Eq. 3 in case of hydration or leaching, and creep computed with Eq. 1 is based on the elastic strain which depends also on Eq. 2. The combination between the different strains is a classical cumulative form:

t T t t t t t i ij i pl ij M ij K ij E ij ij ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ = ∂ ∂ ∑ α δ ε ε ε ε ε , ( 7
)
With superscript E for the elasticity level (Fig. 1), K for the Kelvin level, M for the Consolidating Maxwell level, α the coefficient of thermal expansion and

i pl ij ,
ε the permanent strains associated to each plastic criterion.

Probabilistic scale effect

Experiments on large reinforced structured carried out during the CEO.fr project [START_REF]French National research project for design and assessment of special concrete structure toward cracking and shrinkage[END_REF][START_REF]Comportement et Evaluation des Ouvrages Spéciaux Fissuration-Retrait : Recommandations pour la maîtrise des phénomènes de fissuration ed P des ponts et Chaussées (IREX)[END_REF] have clearly shown the important effect of the concrete heterogeneity on the cracking pattern. In a large structure, the first cracks appear under stresses lesser than the tensile strength measured with usual laboratory tests (such as split test or direct tension test). This is due to the heterogeneity character of the concrete at meso-scale (scale of aggregate), which leads to a random field of tensile strength, and consequently to the existence of weakest links of which the strength decreases when the size of the solicited volume increases. This effect is illustrated in Fig. 5.

Figure 5: Scale effect on tensile strength Rt, versus characteristic size of the specimen (D) and compressive strength fc, adapted from [START_REF] Sellier | Weakest link and localisation WL 2 : a method to conciliate probabilistic and energetic scale effects in numerical models[END_REF] and [START_REF] Rossi | Scale effect on concrete in tension[END_REF].

To consider this effect, the Weibull formulation can be used (Eq. 8); it needs the assessment of the loaded volume V eq which corresponds to the volume of an "equivalent" structure loaded homogeneously; more V eq , lesser R t . In Eq.8, R t,ref is the tensile strength measured with a specimen of which the loaded volume is V ref .

m eq ref ref t t V V R R / 1 ,         = (8) Rt(MPa) fc(MPa) D(mm)
Following the research project CEOS.fr, a new method has been proposed in [START_REF] Sellier | Weakest link and localisation WL 2 : a method to conciliate probabilistic and energetic scale effects in numerical models[END_REF] to actualize the tensile strength at each step of a nonlinear analysis. The proposed method, called WL2, allows considering the strength field heterogeneity without resorting to a random field sampling. Thanks to the WL2 method, the tensile strength used in plastic and damage formulation (Eq.9) takes into account automatically the size of the volume loaded. Indeed, the tensile damage is not assessed with R t,ref , but with R t given by Eq.8.

ij ijkl t ij D σ σ ) 1 ( - = (9)
Informations concerning the damage law are available in [START_REF] Sellier | Orthotropic Damage Coupled with Localised Crack Reclosure Processing[END_REF][START_REF] Sellier | Orthotropic Damage Coupled with Localized Crack Reclosure Processing[END_REF].

Conclusion

The safety demonstration corresponding to the LOCA needs to consider all the phenomena causing either a loss of pre-stressing or a reduction of tensile strength. The best way to merge them is to use a numerical model in which all the phenomena are coupled. The method adopted to clarify the LMDC concrete model consists to build it in the framework of poro-mechanics. The swelling reaction, vapor pressure and water pressure are then considered as internal loadings acting on the solid skeleton. The probabilistic character of the tensile strength is considered through a Weibull formulation (WL2 method); the creep, the plastic strain and the damage depend then automatically on this reduced tensile strength and on the pore pressure. The elastic and the creep strains of the behavior law are written in a differential form. The differential form of the elastics part can be full, or not, depending if the concrete hydrates or dehydrates. The creep velocity derives from the elastic strain to consider implicitly these last chemical effects. The creep strains, corresponding to the nonlinear Maxwell module of the rheological scheme, are obtained thanks to a modification of the characteristic time of creep, which is affected by an anisotropic consolidation coefficient. Once the global architecture of the model implemented, the functions are fitted on several experimental campaigns designed to assess effects of temperature and humidity on the evolution laws associated to this model.

Figure 1 :

 1 Figure 1: Phenomena to consider for a LOCA analysis

Figure 2 :

 2 Figure 2: Poro-mechanic framework to model drying and creep coupling [3] A differential creep model has been proposed by the LMDC [3] to consider both the multi-axiality and the logarithm form of the permanent creep ( M ij ε in Eq.1) in a single differential formulation

)

Figure 3 :

 3 Figure 3: Effects of RH variations between 50% and 98% on the delayed strains of a concrete, f cm = 73MPa, adapted from [4], SL=Simultaneously Loaded and dried, UL=UnLoaded, VD=Variable conditions of Drying, CS=Constantly Sealed Concerning the temperature, the experimental program carried out by Ladaoui [5-7] shows clearly that creep amplitude and creep velocity are both amplified when the temperature rises. This aspect, illustrated in Fig.4 is perhaps among the most important to consider for the LOCA conditions. It is considered in the LMDC model through a dependence of the viscosity (i.e. characteristic time in Eq.1) to the temperature using an Arrhenius law, and through a modification of the creep potential which control the consolidation functions ( ) ( M ij c C ε

Figure 4 :

 4 Figure 4: Influence of temperature on creep velocity and amplitude [5], for two types of HPC with CEMI or CEMV binders Currently, complementary tests are performed by the "Centre de Recherche de l'Industrie du Béton" (CERIB) in the framework of the MACENA research project, in order to quantify creep at higher temperature (until 180°C), and the LMDC carries out tests in various RH and temperature conditions ( 20% ≤ RH ≤ 98%, 20°C ≤ T ≤ 70°C).

Acknowledgement

The authors wish to acknowledge the ANRT for the funding of the project MACENA and EDF for providing the tested materials.


Koen Jochmans
email: koen.jochmans@tse-fr.eu.

"Bias In Instrumental-Variable Estimators Of Fixed-Effect Models For Count Data"

Keywords: count data, bias, fixed effects, inconsistency, instrumental variable, multiplicative-error model, poisson. JEL classification: C23, C26

.

Introduction

The pseudo-poisson maximum-likelihood estimator is routinely used for the purpose of analysing non-negative outcomes. It is consistent under a conditional-mean specification only (Gouriéroux, Monfort and Trognon, 1984a,b) and is well-known to possess a certain robustness against the inclusion of fixed effects (see [START_REF] Wooldridge | Distribution-free estimation of some nonlinear panel data models[END_REF][START_REF] Fernández-Val | Individual and time effects in nonlinear panel data models with large N , T[END_REF][START_REF] Fernández-Val | Individual and time effects in nonlinear panel data models with large N , T[END_REF]. This may lead to the presumption that the instrumental-variable generalizations of the pseudo-poisson estimator, too, are unaffected by the presence of such incidental parameters; it has been used by [START_REF] Tenreyro | On the trade impact of nominal exchange rate volatility[END_REF] and [START_REF] Haucap | How mergers affect innovation: Theory and evidence[END_REF], for example, without any reference to the incidental-parameter problem. This is, however, not the case. This note provides details that underlie this conclusion for the (one-way) panel data case.

Fixed-effect versions of two instrumental-variable estimators are looked at. The first such estimator is the one proposed by [START_REF] Mullahy | Instrumental-variable estimation of count data models: Applications to models of cigarette smoking behavior[END_REF]. The second is the one proposed by [START_REF] Windmeijer | Endogeneity in count data models: An application to demand for health care[END_REF]. These estimators are based on different orthogonality conditions. Furthermore, these conditions are not compatible with one another, a point already made by [START_REF] Windmeijer | Endogeneity in count data models: An application to demand for health care[END_REF]. We show that both estimators are inconsistent under classical asymptotics that treat the length of the panel as fixed, and asymptotically biased under sequences where both dimensions of the panel grow at the same rate. This bias can be corrected for, either by relying on analytical formulae given below or by the application of a jackknife.

The moment conditions underlying the estimator of [START_REF] Mullahy | Instrumental-variable estimation of count data models: Applications to models of cigarette smoking behavior[END_REF] can be modified to yield an alternative estimator that is consistent in short panels. The same would not appear to be true for the poisson type fixed-effect estimator of [START_REF] Windmeijer | Endogeneity in count data models: An application to demand for health care[END_REF].

Panel model and fixed-effect estimators

We observe a scalar outcome y it , a regressor vector x it , and a vector of instruments z it for a random sample of N individuals, i = 1, . . . , N , that are followed over T periods of time, t = 1, . . . , T . Our ambition is to estimate the parameter θ in a multiplicative model of the form

y it = α i λ(x it , θ) v it ,
where λ is a known function, the α i are latent variables that capture any heterogeneity across individuals that does not vary over time, and v it is a time-varying disturbance term.

2

We will treat the α i as fixed effects. Hence, all expectations below are to be understood as being conditional on α 1 , . . . , α N . We will write x i := (x i1 , . . . , x iT ) and z i := (z i1 , . . . , z iT ).

Multiplicative error An instrumental-variable estimator based on E(v it |z i ) = 1, as in [START_REF] Mullahy | Instrumental-variable estimation of count data models: Applications to models of cigarette smoking behavior[END_REF], involves the empirical moment condition

N i=1 T t=1 z it y it α i λ(x it , θ) -1 = 0,
and the the empirical moment conditions for the individual effects. The latter follow from

E(v it ) = 1 and are T t=1 y it α i λ(x it , θ) -1 = 0, for each individual.
Additive error The model can be reformulated as

y it = α i λ(x it , θ) + u it , with u it := α i λ(x it , θ) (v it -1
). An instrumental-variable estimator based on the assumption that E(u it |z i ) = 0, as in [START_REF] Windmeijer | Endogeneity in count data models: An application to demand for health care[END_REF], is constructed around the empirical moment condition

N i=1 T t=1 z it (y it -α i λ(x it , θ)) = 0, together with T t=1 (y it -α i λ(x it , θ)) = 0
for each individual. Clearly, this is an instrumental-variable version of the pseudo poisson estimator.

Large-sample behavior

As discussed by [START_REF] Mullahy | Instrumental-variable estimation of count data models: Applications to models of cigarette smoking behavior[END_REF] and [START_REF] Windmeijer | Endogeneity in count data models: An application to demand for health care[END_REF], the conditions

E(v it |z i ) = 1 and E(u it |z i ) = 0 are not compatible, in general.
The exception is when the regressors are exogenous. Hence, the multiplicative and additive specifications yield (different) estimators based on different conditional validity conditions that will usually not both be valid at the same time. In the context of simultaneous equations the former appears more natural (see, e.g., Mullahy 1997 for a discussion). We treat them separately, in turn.

Multiplicative error It is useful to work with a concentrated estimating equation for θ.

For a given value ϑ of the common parameter the corresponding estimator of α i is available

in closed form. It equals αi (ϑ) := 1 T T t=1 y it λ(x it , ϑ)
.

Substituting this back into the moment condition for θ then yields the profiled estimating equation

N i=1 T t=1 z it y it αi (ϑ) λ(x it , ϑ) -1 = 0.
This estimating equation is biased, in general. Indeed, letting vi :

= T -1 T t=1 v it , its expectation at θ is equal to N i=1 T t=1 E z it E v it vi -1 z i .
This will be non-zero, in general, unless E(v it /v i |z i ) = 1. So, the no-bias condition amounts to conditions on the distributions of the v it |z i that guarantee that the expectation of a ratio is equal to the ratio of expectations. One situation in which this will be the case is when

v it |z i is i.i.d. over time. 1
The presence of a non-vanishing bias in the estimating equation for θ as N → ∞ with T held fixed means that the corresponding instrumental-variable estimator is inconsistent for θ under such asymptotics.

1 In this case v it /v i and v ij /v i have the same distribution (and hence the same mean), conditional on z i .

As this holds for all j, we can average to find that

E(v it /v i |z i ) = 1 T T j=1 E(v ij /v i |z i ) = E(v i /v i |z i ) = 1, as claimed.
The bias does decrease as T → ∞. Observe that αi := αi (θ) is (conditionally) unbiased for the individual effect α i . Furthermore, it is asymptotically linear, as T → ∞, with representation αi -

α i = 1 T T t=1 α i (v it -1) + o p (T -1 ).
A second-order Taylor expansion of the profiled moment condition around the true α i and evaluating at θ yields

N i=1 T t=1 z it (v it -1) - N i=1 T t=1 z it v it α i (α i -α i ) + N i=1 T t=1 z it v it α 2 i (α i -α ) 2 + o p (N ).
The first part corresponds to the infeasible estimating equation for θ in which the individual effects are known and has mean zero. The remaining two terms contribute bias to the feasible equation. If we plug-in the linear representation for αi -α i and take expectations the bias is found to be

- N i=1 T t=1 E(z it v it (v it -1)) T + N i=1 T t=1 E(z it v it ) T T t=1 E(v it (v it -1)) T + o(N ),
where we have assumed (conditional) independence of the errors over time for simplicity of exposition.

The bias in the moment condition translates into a bias of order T -1 in the estimator.

Consequently, while consistent as N, T → ∞, θ will be asymptotically-biased if N and T grow at the same rate. This bias can be estimated and removed, either through the use of a plug-in estimator of the bias based on the analytical formula above, or through a jackknife. Doing so will recenter the limit distribution around zero and lead to asymptotically-valid inference. z it (y it -αi (ϑ) λ(x it , ϑ)) = 0.

Additive errors

When the regressors are strictly exogenous, i.e., when E(u it |x i ) = 0 it is easily seen that E(α i |x i ) = α i . From this, then, unbiasedness of the profiled estimating equation follows, leading to consistency under classical asymptotics. More generally, however, unbiasedness would require that

E T t=1 E(u it |x i , z i ) T t=1 λ(x it , θ) z i = 0
and this will not be true, in general.

We can again characterize the leading bias in the estimating equation. Afirst-order expansion around the α i now suffices as, here, the estimating equation is linear in the individual effects. On the other hand, the estimator of the individual effects is now biased, and this bias has to be accounted for. As T → ∞, we have, again using (conditional) independence over time, αi -

α i = - 1 T 2 T t=1 E(λ it u it ) E(λ i ) 2 + 1 T T t=1 u it E(λ i ) + o p (T -1 ),
where λ it := λ(x it , θ) and λ i := T -1 T t=1 λ it ; this follows from standard higher-order asymptotics [START_REF] Bao | The second-order bias and mean squared error of estimators in time-series models[END_REF]. Together with the expansion

N i=1 T t=1 z it (y it -αi λ it ) = N i=1 T t=1 z it u it - N i=1 T t=1 z it λ it (α i -α i )
this yields the bias as

N i=1 T t=1 E(z it λ it )/E(λ i ) T T t=1 E(u it λ it )/E(λ i ) T - N i=1 T t=1 E(u it z it λ it )/E(λ i ) T + o(N ).
Proceeding as before yields the same conclusions about the large-sample behavior of the instrumental-variable estimator as in the multiplicative case. A bias correction to the estimator may again be constructed.

An alternative estimator

In the multiplicative specification

E(v it |z i ) = 1 implies that E(y it /λ it |z i ) = α i . Therefore, E y it λ(x it , θ) - y it-1 λ(x it-1 , θ) z i = 0
Finding a simple data generating process for the outcome that satisfies the restrictions of the additive specification of [START_REF] Windmeijer | Endogeneity in count data models: An application to demand for health care[END_REF] is not a simple task.

We therefore proceeded as follows. We drew multiplicative errors v it from log-normal distributions depending on the values of the regressor and instrument. To satisfy the condition E(v it |z i ) = 1 we

E(v it |x it , z it ) =              4.4817 if x it = 0 and z it = 0 0.1296 if x it = 1 and z it = 0 0.2231 if x it = 0 and z it = 1 2.8127 if x it = 1 and z it = 1
in a first set of experiments. This is the multiplicative model. To satisfy E(u it |z i ) = 0, in turn, we set

E(v it |x it , z it ) =              4.4817 if x it = 0 and z it = 0 0.6798 if x it = 1 and z it = 0 0.2231 if x it = 0 and z it = 1 1.6669 if x it = 1 and z it = 1
in a second set of experiments. This is the additive model. Tables 1 and2 contain results

for the respective designs for different samples sizes.

Table 1 reports the mean, the standard deviation, as well as the coverage rate of 95% confidence intervals (as obtained over 5000 Monte Carlo replications) for four estimators.

The first, EXO, is the fixed-effect estimator of [START_REF] Mullahy | Instrumental-variable estimation of count data models: Applications to models of cigarette smoking behavior[END_REF] with x it instrumenting for itself. The second, IV, use z it as an instrument for x it . The third, BC, is the (analytically)

bias-corrected version of IV using the formula derived above. (The split-panel jackknife of [START_REF] Dhaene | Split-panel jackknife estimation of fixed-effect models[END_REF] gave very similar results.) Finally, the fourth, DIFF, is the differencing estimator described in the previous section. 2 In all cases coverage rates were computed using the nonparametric percentile bootstrap (based on 99 bootstrap samples), resampling the cross-sectional units. Doing so is particularly useful for the fixed-effect 2 As the covariate is non-negative in our example, the empirical moment condition as stated approaches zero as θ → ∞. Evaluating them at x it -x, where x is the overall sample mean of the x it , in stead of at

x it resolves this problem. Also see [START_REF] Windmeijer | Moment conditions for fixed effects count data models with endogenous regressors[END_REF].

estimators, for which we observed that a plug-in estimator of the asymptotic variance tends to be inaccurate for relatively small values of T . This is in line with observations made elsewhere (see, e.g., Jochmans 2017).

Table 2 has the same structure as Table 1 but concerns the additive specification. Hence, here, all of EXO, IV, and BC are in reference to the [START_REF] Windmeijer | Endogeneity in count data models: An application to demand for health care[END_REF] moment conditions. Results for a DIFF estimator are not reported as such an approach is not available here.

The results show that EXO is heavily biased and not useful for inference in any of the cases considered. This is, of course, expected as it is based on invalid moment conditions.

The IV estimator is biased but consistent as N, T → ∞. This is borne out in our simulations as the bias and standard deviation both shrink in larger samples. However, the bias is important relative to the standard deviation and so hypothesis tests will be size distorted.

In our particular design, this is particularly visible here in Table 1 and less so in Table 2, although it would start becoming more visible also there as the cross-sectional dimension would increase. The corrected estimator, BC, removes most of the bias from IV, re-centering its sampling distribution. For the multiplicative specification in Table 1 we also observe that DIFF is well behaved for all configurations.

Conclusion

This note has highlighted difficulties with instrumental-variable estimators for count data in the presence of fixed effects. The problem can be rectified in long panels (i.e., under rectangular-array asymptotics) and we have shown how to do so. A differencing strategy for short panel data has also been proposed.

Table 1 :

 1 Simulations for multiplicative-error model : instrumental-variable estimator of[START_REF] Mullahy | Instrumental-variable estimation of count data models: Applications to models of cigarette smoking behavior[END_REF] instrumenting x it by itself. IV: instrumental-variable estimator of[START_REF] Mullahy | Instrumental-variable estimation of count data models: Applications to models of cigarette smoking behavior[END_REF] instrumenting x it by z it . BC: bias-corrected version of IV. DIFF: differencing estimator from Section 3. True parameter value: θ = 1.Coverage rates computed through the bootstrap with 99 replications. Results based on 5000 simulations.

				MEAN			STD				COVERAGE
	N	T	EXO	IV	BC	DIFF	EXO	IV	BC	DIFF EXO	IV	BC	DIFF
	200 20 -0.5242 1.3827 1.0683 1.0063 0.0817 0.1893 0.1407 0.1823	0	0.4520 0.9150 0.9440
	200 40 -0.4591 1.1484 1.0249 0.9998 0.0598 0.1028 0.0946 0.1323	0	0.6650 0.9400 0.9330
	200 60 -0.4240 1.0985 1.0165 1.0032 0.0507 0.0791 0.0763 0.1106	0	0.7260 0.9330 0.9340
	200 80 -0.4092 1.0735 1.0107 1.0023 0.0436 0.0655 0.0651 0.0933	0	0.7670 0.9180 0.9230
	200 100 -0.3995 1.0558 1.0051 1.0009 0.0386 0.0586 0.0580 0.0846	0	0.8130 0.9320 0.9400

EXO

Table 2 :

 2 Simulations for additive-error model : instrumental-variable estimator of[START_REF] Windmeijer | Endogeneity in count data models: An application to demand for health care[END_REF] instrumenting x it by itself. IV: instrumental-variable estimator of[START_REF] Windmeijer | Endogeneity in count data models: An application to demand for health care[END_REF] instrumenting x it by z it . BC: bias-corrected version of IV. True parameter value: θ = 1. Coverage rates computed through the bootstrap with 99 replications. Results based on 5000 simulations.

				MEAN			STD			COV
	N	T	EXO	IV	BC	EXO	IV	BC	EXO	IV	BC
	200 20 0.1403 1.0803 1.0049 0.0705 0.2187 0.2038	0	0.9040 0.9250
	200 40 0.1399 1.0456 1.0090 0.0467 0.1367 0.1322	0	0.9150 0.9320
	200 60 0.1421 1.0267 1.0027 0.0379 0.1201 0.1176	0	0.9300 0.9360
	200 80 0.1427 1.0215 1.0037 0.0330 0.1007 0.0991	0	0.9320 0.9390
	200 100 0.1402 1.0018 0.9987 0.0295 0.0890 0.0895	0	0.9320 0.9400

EXO

Support from the European Research Council through grant n o 715787 (MiMo), and from the French Government and the ANR under the Investissements d' Avenir program, grant ANR-17-EURE-0010 is gratefully acknowledged.

estimation from short panels. An example of such an estimator would be the solution to the estimating equation

which is in line with the moment conditions used previously. More generally, optimal (unconditional) moment conditions can be constructed in the usual way [START_REF] Chamberlain | Asymptotic efficiency in estimation with conditional moment restrictions[END_REF]. Under regularity conditions the implied estimator will be N -1/2 -consistent and asymptotically normal as N → ∞ with T fixed. This approach can equally be used in a setting characterized by sequential moment conditions on the form E(v it |z i1 , . . . , z it ) = 1.

In this case

follows from iterating expectations. In both cases, our formulation allows for unrestricted serial dependence in the errors.

For the additive specification where E(u it |z i ) = 0 a differencing strategy would not appear to be available.

Simulations

We generated outcomes using an exponential link function without individual effects, i.e.,

setting θ = 1. The regressor and the instrument were both binary random variables with P(x it = 1|z it = 0) = .80, P(x it = 1|z it = 1) = .30, and P(z it = 1) = .30.


Emmanuelle Petelet-Giraud

Philippe Négrel

Catherine Guerrot

Boron isotope variation during flood events in a Mediterranean basin: Tracer of the water compartments (Hérault, S. France)

Keywords: Boron isotopes, flood event, Mediterranean catchment

Boron isotopes were measured during 3 flood events at the outlet of the Hérault watershed in complement with Sr isotopes. The main water sources (i.e. karst aquifer and the tributaries) were also constrained in high and low flows. Boron isotopes are controlled both by natural and anthropogenic inputs (sewage effluents and/or fertilizers) whereas Sr isotopes trace the various lithologies in the basin. The signature of the Hérault River during flood events presents isotopic variations as large as 15 delta units for boron isotopes, revealing the various contributions of the water compartments according to the rainfall intensity, duration and location.

Introduction

The Hérault River is located in the south of France and drains a medium-sized basin of 150 km long and 2500 km². The watershed is divided into three geological compartments: (1) the upper part is mainly composed by granitic rocks, schist, metalimestone and metadolomite (Palaeozoic basement); (2) in the middle part, Triassic clay and sandstone with some evaporitic layers (anhydrite and gypsum) constitute the Mesozoic cover. The upper layers are formed by Jurassic limestone with some dolomite, clay and calcareous marl. The Jurassic and Cretaceous limestones are folded, faulted and highly karstified; (3) down to the Mediterranean Sea, the alluvial plain is composed of tertiary and quaternary deposits (Fig. 1). The head of the basin is covered by oak, beech and fir trees. In the middle part of the watershed typical Mediterranean vegetation prevails, while most of the alluvial plain is covered with vineyards. The climate is typically Mediterranean and characterized by a highly irregular rainfall regime. Rainfall events are often short but very intense; causing flash floods that can be disastrous. An automatic water sampler device, located at the watershed's outlet in the Mediterranean Sea, collected samples when the river discharge reaches 35 m 3 /s. Between October 2001 and December 2002 and then in January 2006, 7 flood events were automatically sampled with a maximum water discharge between 120 and 1300 m 3 /s, the latter corresponding to the 10-year return flood. Three contrasted flood events were selected for both Sr and B isotopes investigation (Fig. 1).

This study presents δ 11 B results in complement with 87 Sr/ 86 Sr during flood events to decipher anthropogenic versus natural inputs and to assess the anthropogenic impact, mainly from wastewater effluents and/or fertilizers.

Boron concentrations and end-members characterization

The variations in dissolved species content in rivers are particularly noticeable during flood periods with a sharp decrease in concentration with increasing discharge. However, this classical dilution scheme can be highly disrupted in such basin extremely dependent on the flood events and, in a different way, on the different dissolved species considered. In the Hérault watershed, two types of behavior of dissolved species versus the river discharge were evidenced 1 : (1) conservative species like Cl decrease rapidly with increasing discharge according to a fitted logarithmic law; (2) solubility-controlled species like Ca and HCO 3 present no variation with increasing discharge. Boron behaves as a conservative element, as observed in the Seine River 2 , and presents a clear decreasing content with increasing discharge according to a logarithmic law (Fig. 1). It is worth noting that the fitted curve plots above the theoretical dilution line (i.e. dilution with a "zero concentration" water) and defines a plateau with a higher B content than the local rain. As demonstrated for other species 1 (Na, Mg, Sr…), this reflects another B source in addition to the simple rain water dilution. Basically in a river basin, the sources of boron mainly originate from (1) rainwater, (2) water-rock interaction and (3) anthropogenic inputs. The Hérault watershed is a coastal basin mainly recharged by marine rain, i.e. with a marine δ 11 B signature, even if some rain event can have a continental origin with a lower δ 11 B signature. The basin is characterized by various rock types: granite, carbonate, evaporite, clays and volcanic rocks, leading to a potential large range of δ 11 B signature in the draining waters. The anthropogenic activities releasing boron mainly consist in treated urban effluents as boron is widely used in detergent and is not removed by classical waste water treatment and agriculture with vineyards representing about 70% of the cultivated areas [START_REF] Rabiet | Consequence of treated water recycling as regards pharmaceuticals and drugs in surface and groundwaters of a medium-sized Mediterranean catchment[END_REF] and requiring 150-200 g/ha/y of boron for an optimal growing, mostly available in the chemical fertilizers.

At the catchment scale, boron isotopic compositions vary from 5 to 35 ‰ in surface and groundwater, while boron contents vary from 4.7 to 127 μg.L -1 (Fig. 2). The local rainwater signature (δ 11 B=40.3‰; [B]=3.9μg.L -1 ) is typical of coastal rainwater [START_REF] Millot | Multi-isotopic composition (δ 7 Li-δ 11 B-δD-δ 18 O) of rainwaters in France: origin and spatio-temporal characterization[END_REF] . The main karstic spring (H102) and the Vis river (H16) draining carbonates present a slight evolution from rain input both for B content and isotopic signature. From the upstream of the Hérault River down to the outlet of the karst area, the B content is lower than 10 μg.L -1 with a δ 11 B between 22 and 25‰. The tributaries of the alluvial plain present lower δ 11 B together with increasing B contents, that can be related either to anthropogenic inputs with sewage effluents and/or fertilizers, or to natural input, i.e. drainage of argillaceous sediments and evaporites. Three waste water treatment plant effluents were sampled in the alluvial plain which concentrates most of the population. These samples present homogeneous δ 11 B~0‰ typical of sewage effluents with [B]=320±60μg.L -1 .

The Lergue river presents the lowest δ 11 B (5.5‰) and [B]=76μg.L -1 in low flow period (0.95 m 3 .s -1 ), i.e. a boron flux equal to 72mg.s -1 . Two main WWTPs are present in the Lergue basin, S14 releases about 6 mg.s -1 of boron, considering a mean entering discharge 5 of 2000 m 3 /day assumed to represent the output discharge with a [B]=379μg.L -1 , in the same way the second WWTP releases about 10 mg.s -1 of boron. Thus in low flow period (<1m 3 .s -1 ), assuming that all the WWTP effluents reach the river, the B from WTTP can represent a maximum of 20% of the boron flux of the Lergue river. Note that the Lergue River also drains Triassic evaporites that can release B in solution with a δ 11 B > 20‰. The Thongue River presents the highest B concentration with 127 μg.L -1 and δ 11 B = 11.1‰. The Boron flux from the S4 WWTP is about 3 mg.s -1 but cannot be compared to that of the river that is not gauged. Nevertheless, it was demonstrated that this river is impacted by WWTP effluents using the Gd anomaly [START_REF] Rabiet | Consequence of treated water recycling as regards pharmaceuticals and drugs in surface and groundwaters of a medium-sized Mediterranean catchment[END_REF] . The Peyne River was also evidences to be impacted by WWTP effluents [START_REF] Rabiet | Consequence of treated water recycling as regards pharmaceuticals and drugs in surface and groundwaters of a medium-sized Mediterranean catchment[END_REF] .

Boron release from chemical fertilizers is particularly difficult to assess in the catchment as it depends on plant uptake, delay between application period and rain…, nevertheless such an input cannot be excluded and would lead to lower the δ 11 B signature of the main river as δ 11 B of chemical fertilizers are generally between -5 and +15‰.

δ 11 B evolution during a flood event

The flood event of January 2006 was sampled 10 times, including the flood pic at 1150 m 3 /s (Fig. 1B). Sr isotopes were also measured and compared to previous studied flood events [START_REF] Petelet-Giraud | Geochemical flood deconvolution in a Mediterranean catchment (Hérault, France) by Sr isotopes, major and trace elements[END_REF] . The initial signature of the Hérault River at the outlet rapidly evolves when the discharge increases towards that of the upstream carbonated compartment. At the maximum flow (samples 4, 5 and 6, Fig. 2B) the influence of the karstic springs is maximum and then the Hérault signature moves back towards the signature of the alluvial tributaries Lergue and Peyne when the discharge decreases. This behavior is very similar to that observed for previous flood events [START_REF] Petelet-Giraud | Geochemical flood deconvolution in a Mediterranean catchment (Hérault, France) by Sr isotopes, major and trace elements[END_REF] . Similarly to Sr isotopes, the initial δ 11 B signature of the Hérault evolves towards that of the carbonated compartment, i.e. an increase of δ 11 B and decrease of B concentration (Fig. 2C). At the flood pic, all the dissolved boron seems related to the outlet of the karstic area (H11). When the discharge decreases, the δ 11 B decreases again together with an increase of the B concentration, reflecting the contribution of the alluvial tributaries, dominated by the Lergue River. It is worth noting that during the flood, the δ 11 B signature of the alluvial tributaries, especially the main one (Lergue), have to remain low (<10‰) to explain the measured δ 11 B of the Hérault. This thus implies that the natural δ 11 B of the Lergue must be relatively low as the proportion of boron from WWTP inputs becomes very limited (max 200 mg.s -1 ) during high water stages compared to >10 000 mg.s -1 exported at the outlet of the Hérault River. When coupling B and Sr isotopes for 3 flood events (Fig. 2D), both isotopes systematics discriminate the various water compartments of the basin, i.e. the karstic aquifer and the various tributaries of the alluvial plain. During floods, the anthropogenic contributions (flux from WWTP) are negligible compared to natural inputs related to water-rock interaction.

Fig. 1 .

 1 Fig. 1. (A) Simplified geological map and sampling location of the Hérault watershed; (B) Discharge at the outlet of the Hérault River and sampling points; (C) Boron concentration and δ 11 B versus discharge at the outlet of the basin.

Figure 2 :

 2 Figure 2: (A) δ 11 B vs. 1/B, end-member characteristics and Jan-2006 flood, (B) 87 Sr/ 86 Sr vs. 1/Sr during the Jan-2006 flood, (C) δ11B vs. 87 Sr/ 86 Sr for 3 flood events

Acknowledgements

This study was funded by the BRGM Research Division. J. Cubizolles, J-C. Foucher and M. Brach are warmly thanked for their technical support for the conception, installation and maintenance of the automatic sampling device. The authors also thank J.L. Seidel for providing a cumulative rain sample from the studied site.