S Lahaye

T D Huynh

A Tsilanizara

J C Jaboulay

S Bourganel

Comparisons between a priori Uncertainty Quantification and Calculation/Measurement Discrepancies Applied to the MERCI UO 2 Fuel Rod Decay Heat Experiment

INTRODUCTION

In 2008, an UO 2 fuel rod was irradiated up to ∼ 3.5GWd/t HL in the CEA/Saclay research reactor OSIRIS [START_REF] Loubière | OSIRIS Reactor for Fuel, Material and Other Irradiation[END_REF] through the MERCI experiment [START_REF] Blandin | MERCI-MOSAIC: Experimental Tools for Residual Power Measurement in OSIRIS Reactor[END_REF]. Some rod's pellets were analysed. Experimental results were obtained:

• fuel rod decay heat measurement by calorimetry (cooling times from 27 minutes to 42 days),

• evaluation of the amounts of some nuclei (U, Pu, Cs, Nd...) by isotopic dilution mass spectrometry,

• nuclide activities.

Comparison of type (C -M)/M where C is the calculated decay heat and M the experimentally measured decay heat was reported in reference [START_REF] Jaboulay | Analysis of MERCI Decay Heat Measurement for PWR UO 2 Fuel Rod[END_REF]. Neutronics flux was computed by CEA TRIPOLI-4 Monte-Carlo or APOLLO2 deterministic transport code system [START_REF] Brun | Systematic Uncertainty Due to Statistics in Monte Carlo Burnup Codes: Application to a Simple Benchmark with TRIPOLI-4-D[END_REF][START_REF] Sanchez | APOLLO2 year 2010[END_REF]. Irradiation and cooling phases were simulated with CEA DARWIN/PEPIN2 [START_REF] Tsilanizara | DARWIN: an Evolution Code System for a Large Range of Applications[END_REF] code system. Good agreements were globally shown between numerical simulation results and experimental measurements on total decay heat. Uncertainty propagation in this previous work [START_REF] Jaboulay | Analysis of MERCI Decay Heat Measurement for PWR UO 2 Fuel Rod[END_REF] was established from Rebah's studies [START_REF] Rebah | Uncertainty on Decay Heat Power Due to Fission Product Data Uncertainty[END_REF].

These last years, within the framework of supports to its industrial partners EDF and Areva, the CEA/DEN has implemented two different approaches to propagate nuclear data uncertainties to depletion code outputs. A direct forward deterministic perturbation method using sensitivity profiles was implemented in the French industrial depletion code system DARWIN/PEPIN2 [START_REF] Tsilanizara | DARWIN: an Evolution Code System for a Large Range of Applications[END_REF]. Meanwhile, a propagation based on Monte Carlo correlated sampling [START_REF] Tsilanizara | Probabilistic approach for decay heat uncertainty estimation under URANIE platform by using MENDEL depletion code[END_REF] was achieved using the CEA uncertainty platform URANIE [START_REF] Gaudier | URANIE: the CEA/DEN Uncertainty and Sensitivity platform[END_REF] and the new generation depletion code MENDEL [START_REF] Lahaye | First Verification and Validation Steps of MENDEL Release V1.0 Cycle Code System[END_REF]. Both methods will be described in this article. Uncertainties on nuclear data (independent fission yields, multigroup microscopic cross sections, radioactive decay constants, radioactive decay branching ratio and radioactive decay energies) are propagated to decay heat and isotopic concentrations.

To match MERCI first study [START_REF] Jaboulay | Analysis of MERCI Decay Heat Measurement for PWR UO 2 Fuel Rod[END_REF], flux computation is realized by APOLLO2. Furthermore, uncertainty quantification (UQ) is re-evaluated by both DARWIN/PEPIN2 and MENDEL codes. It is an important stage in the process of validating both depletion codes with their respective data library.

UNCERTAINTIES ON NUCLEAR DATA

Experimental results are given with a discrepancy, which is assimilated to the standard deviation of the random variable associated to the measurement value.

Uncertainty quantification on numerical simulations aims to compute the standard deviation of the outputs due to the uncertainties on the input data. In this paper, we consider only uncertainties due to nuclear data. Other potential sources of uncertainties are considered perfectly known (geometry, technological data...). Nuclear data without uncertainty data in the evaluation files are associated to a zero standard deviation (no uncertainty).

In the present work, DARWIN/PEPIN2 and MENDEL use uncertainty data from JEFF-3.1.1 [START_REF] Kellet | [END_REF] for independent fission yields, radioactive decay constants, radioactive decay branching ratios and radioactive decay energies. Table I shows the quantity of uncertain parameters taken from JEFF-3.1.1. No correlations are considered between different physical quantities. Radioactive decay periods and radioactive decay energies are two by two independent. Radioactive decay branching ratios are correlated in such a way that the sum of branching ways is equal to 1 (correlation matrix).

For independent fission yields, MENDEL and DARWIN/PEPIN2 can choose different models: no correlation between yields, normalization (sum of yields equal to a constant for each fissile systems) or correlations taking into account the mass yields constraints [START_REF] England | ENDF-349, Evaluation and Compilation of Fission Product Yields[END_REF].

Uncertainties for microscopic multigroup cross sections are taken from the COMAC data base [START_REF] Al | Estimation of Multi-Group Cross-Sections Covariances for 235[END_REF]. COMAC matrices describe correlations between partial cross sections for one given isotope, and between groups for one given cross section. Data between two distinct nuclei are not correlated. Multigroup data are directly taken into account in URANIE sampling, and multigroup samples are used in MENDEL, while one-group variance and correlations on monocinetic reaction rates are calculated to be used by DARWIN/PEPIN2. It has been shown through inter-comparisons that the effect of the multigroup mesh rather than one group data is of second order.

UNCERTAINTY PROPAGATION METHODS

We describe briefly the propagation methods and hypotheses in this section. Further details will be given in full presentation.

Deterministic Method Used in DARWIN/PEPIN2

DARWIN/PEPIN2 propagates uncertainties from nuclear data to decay heat or isotopic concentrations using a direct forward first order perturbation method. With X the uncertain input variables and Y the uncertain outputs, we can use the following formula:

cov(Y) = S Y/X cov(X)S T Y/X (1)
where cov(X) (resp. cov(Y)) stands for the variancecovariance matrix for variable X (resp. Y) and S Y/X stands for the sensitivity matrix of Y regarding X.

The use of a direct forward first order perturbation method implies the hypothesis of linearity of the outputs (as a function of uncertain nuclear data). For small perturbations of the uncertain input parameters (limited to one standard deviation), this approximation is proved to be valid when comparing with propagation approaches not taking this linearity hypothesis, such as stochastic uncertainty propagation approach used in MENDEL.

The uncertainty propagation in DARWIN/PEPIN2 needs the same number of depletion calculations as the number of uncertain parameters for the calculation of sensitivity coefficients. This is performed in parallel mode by DARWIN/PEPIN2 through its INCERD module, which establish the sensitivity matrix before computing the interest output covariance through equation ( 1).

Stochastic Method Used in MENDEL

MENDEL propagates uncertainty using a correlated sample method. All realizations are generated with URANIE.

For this study, 2000 realizations generated by LHS sampling method were propagated. Most random variables are supposed to be Gaussian if the relative standard deviation is lower than 50% and are supposed to be Log-Normal in the other cases. Fission yields are always considered log-normal. Correlations are taken into account when available, through the use of importance sampling in URANIE.

The sampling method and hypotheses will be described in the full presentation.

DECAY HEAT UNCERTAINTY QUANTIFICATION

We compare numerical calculation discrepancy to decay heat measurements and uncertainty propagation results (likelihood) in Figs 1 and2. The red plain lines show the C -M M relative discrepancy between mesurement (M) and calculations (C). For those computed values, uncertainty on nuclear data has not been taken into account. Grey domain marks an uncertainty zone of 2 standard deviation around the reference value (zero). Note that the decay heat measurement uncertainty is given to ±1% and is included to the grey domain definition.

This grey domain corresponds to a 97.7% likelihood if decay heat is considered to be Gaussian.

Both graphs are similar, and prove a rather good agreement between experimental data and computation. Indeed, for all considered cooling times (from 27 min to 42 days), the two standard deviation domain obtained by uncertainty quantification covers the decay heat discrepancy between experimental data and calculation results. Contributions by type of nuclear data are shown in Fig. 3. The main contributor to decay heat uncertainty among nuclear data are the uncertainties due to independent fission yields (in red). We considered here for both codes that fission yields are normalized so that the sum for one fissile system is constant. We do not take here in consideration the correlations resulting from mass yields constraints [START_REF] England | ENDF-349, Evaluation and Compilation of Fission Product Yields[END_REF]. For this small burnup, other contributors are negligible, especially radioactive decay branching ratios and radioactive decay periods.

Full presentation will analyze more precisely uncertainty quantification and sensitivity analysis, in particular the effects of correlation hypotheses on fission yields.

ISOTOPIC DENSITY UNCERTAINTY QUANTIFICA-TION

Experimental ratios of concentrations were measured for both heavy nuclides (U, Pu) and fission products (Cs, Nd). Those values, associated to experimental measurement error bars, enable to verify the actual burnup of the fuel rod, in particular when considering Neodymium isotopes.

We show in graph 4 experimental values and MENDEL uncertainty quantification propagation values for Neodymium ratios. DARWIN/PEPIN2 uncertainty propagation results are coherent with MENDEL ones for all isotopes. Error bars correspond to one standard deviation. We observe a good agreement for all ratios, all experimental error bars being fully included in MENDEL uncertainty quantification standard deviation.

For Plutonium, only Pu242/Pu239 experimental value is outside the one standard deviation of calculation results, due to the very low amount of Pu242, leading to a bad evaluation of the quantities during the isotopic dilution mass spectrometry.

For Cesium and Uranium isotopes MERCI experimental results are inside the one standard deviation of calculation results, as will be shown in full presentation.

CONCLUSION

Comparisons between codes and experimental measurement give globally good agreement.

Other isotopic concentration uncertainty propagation re-sults and analysis of correlation hypotheses consequences on decay heat uncertainty will be given in the full presentation. This work is a new contribution for the validation of both MENDEL and DARWIN/PEPIN2 codes for isotopic concentration and decay heat computations, as well as both deterministic and probabilistic methods to propagate nuclear data uncertainty to those quantities.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig.1. In red, decay heat discrepancy between experiment data and DARWIN. In grey, the 2 standard deviation from deterministic method used in DARWIN/PEPIN2 domain.

Fig. 4 .

 4 Fig. 4. MERCI experimental data and MENDEL uncertainty quantification for Neodymium isotopes.

TABLE I .

 I Number of physical parameters with non-zero uncertainty value in JEFF-3.1.1 evaluation.

	evaluation	JEFF-3.1.1
	Independent fission yields	33668
	Radioactive decay constant	3204
	Radioactive branching ratios	505
	Radioactive decay energies	1554

ACKNOWLEDGMENTS

The authors thank EDF and AREVA for their financial support to the developments of DARWIN/PEPIN2 and MENDEL code systems.


C Cavadias

G Mariet

Cote

    Predict and compare the formation of segmented flow in microsystems: Interest for radiochemical liquid-liquid extraction

Microscale assets [1]  Volume reduction (from mL to µL)  Lower operators' exposure to  Shorter analysis time  Lower polydispersity  Cost reduction Vansteene A. 1 , Jasmin J.P. 1 , Vera de la Cruz G. 2 , Brennetot R. 1 , Cavadias S.

PREDICT AND COMPARE THE FORMATION OF SEGMENTED FLOW IN MICROSYSTEMS : INTEREST FOR RADIOCHEMICAL LIQUID-LIQUID EXTRACTION

One of the most important separation techniques in radiochemical procedures is solvent extraction. In the last decade, a growing interest in its use in microsystems has emerged because such systems allow a good control of both the interface area between aqueous and organic phases, and the contact time of the two phases. To improve mass transfer in liquid-liquid extraction reactions for radiochemical separation relevant to the nuclear field, our goal is to implement liquid-liquid extraction with segmented flows in microsystems in order to increase the interfacial area. The optimization of the interfacial area requires taking into account the hydrodynamics of biphasic systems. This study is dedicated to the understanding of segmented flows in a T-junction with the aim to form tailor-made droplets for given organic and aqueous phases.

Based on previous works from Van Steijn and Xu, a numerical model was developed to accurately predict the characteristics of droplet production at a T-junction (volumes, frequency, spacing, and specific interfacial area). This numerical model was verified experimentally by generating on-chip segmented flows, and offered ways of improvement in both operating parameters and chip design. Further works will focus on porous hydrophobic membrane (PTFE) phase separation and on implementation of droplet generation for relevant nuclear chemical systems. This work will enable us to characterize each phase during liquid-liquid extraction in order to take into account both matter exchange and hydrodynamics in following mass transfer studies.

MODEL VALIDATION FOR H 2 O/DODECANE CHEMICAL SYSTEM WITH DROPLET MORPHOMETRY AND VELOCIMETRY (DMV) ANALYSIS

Photron AX100 , recording at 6.000 fps, replay at 30 fps Thanks to a video processing software for timeresolved, label-free tracking of droplet parameters [4] , we determined experimentally the characteristics of the generated flow.

P dodecane = 143600 Pa, then Q dodecane = 2,07 mL.h -1 P water = 153000 Pa, then Q water = 7,22 mL.h -1

Droplet volume comparison between the numerical model based on equations from [2] and [3] and experimental results for the water/dodecane chemical system.

ASSETS OF SEGMENTED FLOWS IN MICROSYSTEMS

diffusion (by using droplets with a small spacing in between) convection (by recirculating species inside droplets)

HYDRODYNAMICS OF SEGMENTED FLOWS IN A T-JUNCTION

Hydrodynamic regimes

Capillary number 𝑪𝒂 𝒄 = η 𝒄 𝑸 𝒄 σ𝒘 𝒄 𝒉 with 𝜂 𝑐 : viscosity of the continuous phase

𝜎: interfacial tension between the aqueous and organic phases

To form tailor-made droplets for given organic and aqueous phases, we identified the geometrical parameters of the T-Junction and the physico-chemical parameters of the aqueous and organic phases involved in the control of the flow regimes in addition to the flow rates.

Aqueous phase Organic phase

For the different regimes, we found equations linking flow rates and physico-chemical parameters of the biphasic chemical system with characteristics of the segmented flow (droplet volume, spacing, …). Therefore, we are able to predict the characteristics of the generated flow.

[4] Basu, A. Lab Chip, 2013, 13, 1892-1901.

  2 , Mariet C. 1 , Cote G. 3 axel.vansteene@cea.fr; jean-philippe.jasmin@cea.fr 1) Den -Service d'Etudes Analytiques et de Réactivité des Surfaces (SEARS), CEA, Université Paris-Saclay, F-91191, Gif sur Yvette, France 2) UPMC -Univ Paris 06, 4 Place Jussieu, 75005 Paris, France 3) PSL, Chimie ParisTech -CNRS, Institut de Recherche de Chimie Paris, 75005, Paris, France

  Fluid Interfaces -Lyon, October 12-13, 2016 [1] Helle G et al. Talanta, 2015, 139: 123-131. Organic phase Aqueous phase COMSOL® 2D CFD Model developed by our team showing evidence of species recirculation inside 70:30 (v/v) dodecane/TBP droplets, in a continuous phase of [HNO 3 ] = 3 M /[U(VI)]=10 -4 M, with a pressure gradient of 0.66 Pa.μm -1 .Van Steijn, V. et al. Lab Chip, 2010, 10(19), 2513-2518. Xu, J. H. et al. Microfluidics and Nanofluidics, 2008, 5(6), 711-717.𝑤 𝑑 = 𝑤 𝑐 = 𝑤 𝑜 = 105 𝜇𝑚 ℎ = 100 𝜇𝑚 𝑤 𝑜 ′ = 300 𝜇𝑚



  Droplet volume 𝑉 𝑝𝑙ug or length l 𝑝𝑙ug  Droplet spacing 𝑒 = 𝑄 𝑐 +𝑄 𝑑 ℎ𝑤𝑓 Generation


Elodie Chanthery
email: elodie.chanthery@laas.fr

Audine Subias
email: audine.subias]@laas.fr

Diagnosis approaches for detection and isolation of cyber attacks and faults on a two-tank system

Recently, a two-tank benchmark for detection and isolation of cyber-attacks has been proposed to the diagnosis research community to test different diagnosis methods. In this work, we use this benchmark, add some scenarios, and test several diagnosis techniques to evaluate their diagnosability power. We propose to use a well-known model-based diagnosis approach to identify the groups of faults and attacks that are not isolable or detectable. In a second part, we investigate improvements provided by data-based diagnosis techniques and show that they succeed to isolate faults and attack for this benchmark.

Introduction

Connected systems and intelligent systems are largely widespread in the majority of industrial domains (production, transport, networks, aeronautics, space,...). Critical systems as power plants, water plants, smart grids also rely heavily on these Information and Communication Technologies (ICT) to control, to monitor and to manage remotely several equipments. If ICT has contributed to improve the control and monitoring of installations it has also exposed critical systems to a new kind of threats by introducing a high degree of connectivity. Consequently, these systems which include both computational resources, communication capabilities and hardware also called cyber-physical systems, are exposed to malicious attacks, known as cyberattacks. As in any computer system connected to internet the well known CIA (Confidentiality, Integrity, and Availability) triad of a cyber-physical system can be threatened by non-targeted attacks, but it can also fall prey of targeted attacks. In this case the hacker knows that he is targeting a control system and the cyber-attack is a malicious act aimed at degrading or interrupting the operation of an industrial installation or even destroying it [START_REF] Cardenas | Attacks against process control systems: Risk assessment, detection, and response[END_REF]. These attacks exploit the vulnerabilities of the system and can have a very important impact not only from a process point of view but also from financial and environmental ones. Obviously, the influence the hacker can have on the system depends on his knowledge about the system. A prominent example is the attack on the water distribution system in Queensland, Australia [START_REF] Slay | Lessons learned from the maroochy water breach[END_REF]. The hacker used a laptop and took the control of 150 sewage pumping stations to release one million liters of untreated sewage into a storm-water drain connected to local waterways. Another example is the Stuxnet attack that taking the control of actuators and sensors struck the Iranian nuclear facility causing the failure of the centrifuges. More recently, we can cite the cyber-attack against Ukraine power grid that cut power to many customers for about six hours. Histories and descriptions of cyber-attacks can be found in [START_REF] Mclaughlin | The cybersecurity landscape in industrial control systems[END_REF][4] [START_REF] Miller | A survey scada of and critical infrastructure incidents[END_REF] [START_REF] Khorrami | Cybersecurity for control systems: A process-aware perspective[END_REF].

In the context of a classical closed-loop control system architecture, the problem of attack diagnosis is getting closer to the problem of fault diagnosis. Nevertheless, the attack diagnosis problem addresses a larger spectrum of feared scenarios. Indeed, attacks as faults can concern sensors, actuators or the controller itself but the hacker is endowed with intelligence and has capabilities that multiply the possibilities of attacks. Whatever the attack type, its effect is at least one of the following: the sensor returns incorrect measurements, the controller sends out incorrect signals, the actuator ignores control input or executes incorrect commands. According to [START_REF] Cárdenas | Research challenges for the security of control systems[END_REF] two main types of attacks can be distinguished: those that concern data integrity issued from sensors or controllers (called by the authors deception attacks) and those of type Denial of Service (DoS) that threat the signals availability preventing the controller from receiving sensor measurements or the actuators from receiving control commands. In the first case, signals are corrupted by the hacker so that the measurement is incorrect, delayed, replayed etc.

The growing interest in this problem of attack detection and diagnosis is reflected in the Diagnosis community. In [START_REF] Amin | Cyber security of water scada systems-part i: Analysis and experimentation of stealthy deception attacks[END_REF] an observer based method is developed to diagnose stealthy deception attack that can enable remote water pilfering from automated canal systems. The work of [START_REF] Pasqualetti | Attack detection and identification in cyber-physical systems[END_REF] define the notions of detectability and identifiability of an attack by its effect on output. [START_REF] Amine | Medium access scheduling for input reconstruction under deception attacks[END_REF] proposes an approach to design an input reconstruction filter to identify the effect of random deception attack on the system behavior in the context of limited access to sensor measurements. In [START_REF] Helem | Detection of replay attacks in cyber-physical systems using a frequency-based signature[END_REF] the authors propose a frequency based approach to detect replay attacks on cyber-physical systems. [START_REF] Sicard | Critical States Distance Filter Based Approach for Detection and Blockage of Cyberattacks in Industrial Control Systems[END_REF] propose an approach based on the introduction of a filter between the controller (Programmable Logic Controller) and the sensors/actuators to detect and prevent orders that could damage the system. Additional information issued from a notion of distance between states allow distinguishing cyberattacks from classical failures. [START_REF] Kawakami Carvalho | Detection and mitigation of classes of attacks in supervisory control systems[END_REF] consider the problem of intrusion detection and mitigation in supervisory control systems. Several types of attacks are considered (en-Figure 1: two-tank system abling/disabling of actuators events, erase/insertion of sensor events).

Recently a two-tank system well known in the diagnosis community [START_REF] Bouamama | Diagnosis of a two-tank system[END_REF] has been extended and proposed to the scientific community to test different approaches of detection and isolation of cyber-attacks. To initiate the work the authors considered a classical model-based approach (Analytical Redundancy Relations) [START_REF] Quevedo | A two-tank benchmark for detection and isolation of cyber attacks[END_REF]. This paper is an answer to this invitation. The contribution of the paper is to show that the initial results obtained from Analytical Redundancy Relations can be improved by the use of data-based methods.

This paper is organized as follows. Section 2 presents the two-tank benchmark. Section 3 describes the different existing scenarios of faults and attacks, as well as some new scenarios proposed in this work. Model-based Diagnosis techniques are applied on the two-tank in Section 4 for detecting and isolating faults and attacks. The use of Data-based methods is developed in Section 5. Section 6 concludes the paper.

Benchmark Description

The two-tank benchmark is illustrated in Figure 1. The goal of the system is to provide water to customers with a continuous flow Q o . A valve V o simulates the consumer. V o is opened in nominal mode.

The system consists in two tanks T 1 and T 2 in series. h 1 and h 2 are supposed to be the water level in T 1 and T 2 , respectively. T 1 is filled by a pump P 1 controlled by a PI level controller acting on the inlet flow Q p . It is assumed that Q p is proportional to the PI controller output U m p , that is measured. Note that in the article, the letter "m" is added to a variable when it is measured. The water flow between the two tanks is Q 12 . It can be controlled by an ON/OFF valve V b . The water levels in the two tanks may be affected by real leaks but also by attacks simulating a leak. V f1 (resp V f2 ) denotes a leakage in T 1 , (resp. in T 2 ) while Q f1 (resp. Q f2 ) is an attack. Two sensors provide the measurements of the real water levels of the tanks: h m 1 (for water level in T 1 ) and h m 2 (for water level in T 2 ). A detailed description of the dynamic model of the faultless two-tank system is given in [START_REF] Bouamama | Diagnosis of a two-tank system[END_REF].

The input u of the system is supposed to be known.

u =    Q m p U m b U m 0 U m p    (1)
where Q m p is the measured flow from the pump

P 1 , U m b represents the position of valve V m b (0 if the valve is open, 1 if the valve is closed), U m
o is the position of the valve V o , and U m p corresponds to the PI controller output.

Fault and attack scenarios

The next sections presents the faults and attack scenarios implemented in the two-tank system benchmark.

Fault Scenarios

A set of fault scenarios has been described in [START_REF] Bouamama | Diagnosis of a two-tank system[END_REF]. f 0 faultless mode: the process runs without fault.

f 1 Pump fault: from 40s up to 120s the pump is simulated off (like if it is broken). f 2 Level sensor mh 1 stuck at zero fault: from 40s up to 120s the output of the sensor is stuck at zero. f 3 Level sensor mh 2 fault: from 40s up to 120s the output of the sensor is stuck at zero. f 4 Leakage fault in Tank T 1 from 40s up to 120s.

Q f1 = 10 -4 m 3 /s. f 5 Leakage fault in Tank T 2 from 40s up to 120s. Q f2 = 10 -4 m 3 /s.
f 6 Sensor fault for the pump P 1 flow: the sensor is stuck at zero from 40s to up to 120s.

f 7 Valve V b stuck-closed fault: v b is stuck closed from 40s
to up to 150s. f 8 Valve fault: U m b is stuck at zero from 40s to 120s. Some additional fault scenarios have been defined and tested in this article.

f 9 Pump sensor fault: Q m p is stuck at zero from 40s to up to 120s. f 10 PI controller fault: U p is always equal to zero from 40s to up to 120s.

f 11 Valve V b stuck-open fault: v b is stuck open from 40s to
up to 150s. f 12 Level sensor mh 1 stuck at non-zero fault: from 40s up to 120s the output of the sensor is stuck at a non-zero value, here h m 1 = 0.4.

Attack Scenarios

In this article we consider as in [START_REF] Quevedo | A two-tank benchmark for detection and isolation of cyber attacks[END_REF] that the final objective of an attacker is to steal water but we consider also the case where the aim of the attacker is the destruction of the installation or the degradation of the behaviour. We consider an active attacker who has the ability to altering the cyberphysical system by exploiting system security design holes.

The attacks can then have some effects on the sensors so that they return incorrect measurements, on the actuators to execute incorrect control actions but also on the controller to perform a physical action on the two-tank system. We recall here the set of attack scenarios defined in [START_REF] Quevedo | A two-tank benchmark for detection and isolation of cyber attacks[END_REF].

a 1 Short-term water theft from T 1 : this scenario is similar to the scenario with fault f 4 . The difference is that it is cast maliciously, with the purpose of stealing water from T 1 . V f1 = 10 -4 m 3 /s from 40s up to 80s. a 2 Short-term water theft from T 1 with hiding signal added to the measurement from 40s up to 80s. The theft is hidden by adding a signal to the output of the level sensor in tank T 1 . The PI controller works as if nothing had happened. a 3 Long-term water theft from T 1 with hiding signal added to the measurement: the attack is the same as a 2 but its duration is extended so that Tank T 1 might become empty, affecting the tank T 2 and Q 0 due to interconnection. a 4 Long-term water theft from T 1 with small signal added to the measurement: the attack is the same as a 3 but the added signal only compensates 50% of the stolen water, so that the pump will compensate half the theft and the attack will be harder to detect. a 5 Short-term water theft from T 2 (same as a 1 ) V f2 = 10 -4 m 3 /s from 40s up to 80s. a 6 Short-term water theft from T 2 with hiding signal added to the measurement from 40s up to 80s (see a 2 ). a 7 Long-term water theft from T 2 with hiding signal added to the measurement (see a 3 ). a 8 Long-term water theft from T 2 with small signal added to the measurement (see a 4 ). a 9 Replay attack from 160s to up to 200s: the hacker records the measurements coming from the sensors without stealing water from the tanks. Then, when the system has reached its steady-state, the hacker steals water while replacing the sensor values by the recorded ones. A set of five new attacks has been added to this initial set as follows: a 10 Overflow with sensor hack from 40s up to 120s: the objective is not to steal water but to destroy the system. The value of h m 1 is modified to simulate a fake leakage, so that the pump compensates and causes a water overflow in T 1 . a 11 Overflow with controller hack from 40s up to 120s: the hacker disconnects the PI controller and sets a too high water flow so that there is an overflow in T 1 . a 12 Sensor hack without theft in tank T 2 from 40s up to 120s: the objective is not to steal water but disturb the system. The sensor value in T 2 is corrupted so that h m 2 remains constant, while the tank is emptied by the customers. a 13 Water theft from T 2 when V 0 is closed from 40s to up to 120s: when Q 0 = 0, the hacker pretends to be a user but the flows are different. a 14 Water theft from T 2 when V 0 is closed from 40s to up to 120s: when Q 0 = 0, the hacker pretends to be a user and the flows are the same.

Model-based diagnosis techniques on the benchmark

In this article we extend the diagnosis results obtained on the two-tank system based an Analytical Redundancy Re-lations [START_REF] Staroswiecki | Analytical redundancy relations for fault detection and isolation in algebraic dynamic systems[END_REF]. Nevertheless, it would be interesting to investigate other model-based approaches notably to take advantage of the dynamical information embedded in the fault and attack profiles. By considering the hybrid nature of the twotank system, model based approaches developed in the field of switched systems could be investigated as observer based methods.

Analytical Redundancy Relations

A set of four Analytical Redundancy Relations (ARR) has already been described in [START_REF] Bouamama | Diagnosis of a two-tank system[END_REF] and four residuals have been included in the benchmark. We aim to use them for fault detection and isolation, but also for attack detection and isolation.

As a reminder, the ARR are calculated as:

r 1 (t) = -C vb U m b (t)sign(h m 1 (t) -h m 2 (t)) |h m 1 (t) -h m 2 (t)| +Q m p (t) -A 1 dh m 1 dt r 2 (t) = C vb U m b (t)sign(h m 1 (t) -h m 2 (t)) |h m 1 (t) -h m 2 (t) -Cvo h m 2 (t)U m 0 (t) -A 2 dh m 2 dt r 3 (t) = U m p (t) -K P (h 1,ref -h m 1 (t)) -K I (h 1,ref -h m 1 (τ ))dτ r 4 (t) = Q m P (t) -          U m p (t) if 0 < U m p (t) < qp,max 0 if U m p (t) ≤ 0 Qp max if U m p (t) ≥ qp,max (2)
with: C vb the hydraulic flow coefficient of the valve V b , A i,(i=1,2) the cross-section of the cylindric tank T i , C vo the hydraulic flow coefficient of the valve V o , K P and K I the coefficients of the PI controller. The numerical values of these parameters are given in [START_REF] Quevedo | A two-tank benchmark for detection and isolation of cyber attacks[END_REF] and [START_REF] Bouamama | Diagnosis of a two-tank system[END_REF].

Two additional residuals have been included in the benchmark, supposing that the pressures in tank T 1 and T 2 may be measured by two new sensors.

We then defined:

r5(t) = P m 1 (t) -ρ.g.h m 1 (t) -Patm r6(t) = P m 2 (t) -ρ.g.h m 2 (t) -Patm (3)

Definition of detectability and isolabity for faults and attacks

Let F = {f 1 , . . . , f n f } be the set of faults of the system, A = {a 1 , . . . , a n a } be the set of attacks. Let f 0 corresponds to the nominal mode without any fault or attack. We define in this section some fundamental definitions for fault diagnosis and attack detection and isolation.

We recall here the definition of a fault signature matrix and extend it for the attacks. Definition 1 (Fault Signature). Given a set ARR composed of n r ARR and F the set of considered n f faults for the system; consider the function ARR × F -→ {0, 1}, then the signature of a fault f ∈ F is the binary vector F S(f ) = [τ 1 , . . . , τ n r ] T where τ k = 1 if f is a variable of the equation used to form arr k ∈ ARR, otherwise τ k = 0.

The same definition holds for attacks. Definition 2 (Attack Signature). Given a set ARR composed of n r ARR and A the set of considered n a attacks for the system; consider the function ARR × A -→ {0, 1}, then the signature of an attack a ∈ A is the binary vector F S(a) = [τ 1 , . . . , τ n r ] T where τ k = 1 if a is a variable of the equation used to form arr k ∈ ARR, otherwise τ k = 0.

We can then define the fault and attack signature matrix as follows: Definition 3 (Fault and attack signature matrix F ASM ). The signatures of all the faults in F and of all the attacks in A together constitute the fault and attack signature matrix F ASM for the system, i.e F ASM = [F S(f 1 ), . . . , F S(f n f ), F S(a 1 ), . . . , F S(a n a )] T . This matrix obviously implies what is usually called the fault signature matrix and what could be called the attack signature matrix.

From the F ASM detectability and isolability properties can be defined for faults and attacks. Definition 4 (Detectable fault). A fault f ∈ F is detectable in the system if F S(f ) = F S(f 0 ).

The same definition can be applied for attacks. Definition 5 (Detectable attack). An attack a ∈ A is detectable in the system if F S(a) = F S(f 0 ). Definition 6 (Isolability). Two faults or attacks f a ∈ F ∪ A and f a ∈ F ∪ A are isolable in the system if they are detectable and if F S(f a) = F S(f a ). Definition 7 (Diagnosability group). Two faults or attacks f a ∈ F ∪ A ∪ {f 0 } and f a ∈ F ∪ A ∪ {f 0 } are in the same diagnosability group if they are not isolable from each other.

In others words, a group of diagnosability is defined as a set of behaviours (nominal, faulty or attacks behaviours) that can not be isolated two by two as their signatures are identical.

Fault and attack signatures

Some diagnosis results have been provided in [START_REF] Bouamama | Diagnosis of a two-tank system[END_REF] but no link between the fault scenarios and the residuals have been clearly established. Here, we aimed at giving the fault and attack signatures for each defined scenario.

Using the four first residuals we succeeded in defining eight diagnosability groups as presented in Table 1.

From Table 1, we can deduce that all faults and attacks are detectable, except the replay attack. Only 3 faults (f 1 ,f 6 and f 10 ) are isolable from the others and from the attacks. Attacks are always included in diagnosability groups with at least one fault so it is impossible to distinguish them from faults with the ARR approach.

We also tried to take benefit from the two additional residuals obtained from pressure sensors (see section 4.1). From Table 2, we can deduce that all faults and attacks are detectable, except the reply attack, as with 4 residuals. The 2 new residuals increase isolability power: for example diagnosability group G 2 has been split into 3 other diagnosability groups G 21 , G 22 and G 23 , while G 4 is divided into 2 diagnosability groups. The only interesting point is that G 42 only refers to attack scenarios a 6 , a 7 , a 8 and a 12 so that it is possible to distinguish fault from attack in these cases, moreover f 3 is now isolable. However, much ambiguity remains between faults and attacks in 5 diagnosability groups (G 21 , G 3 , G 41 G 6 , G 8 ).

Results

An HMI has been provided in the matlab benchmark to visualize online the residuals value and the possible diagnoses. For example, we simulated scenario f 9 and found the results illustrated on Figure 2. r 4 is equal to 0, then rises to 1 at 41s. r 1 is equal to 0, then rises to 1 at 44s. The other residuals remain equal to zero. Consequently, the diagnosis hypotheses are the following: from 0s to 41s, the diagnosis hypothesis in {f 0 , a 9 }. Between 41s and 44s, it is {f 1 }, then from 44s, the final diagnosis hypothesis is {f 9 , a 11 }. It can be noticed that the detection delay is then about 1s and the isolation delay is about 4s in this case.

Data-based Diagnosis techniques on the benchmark

Since model-based methods are not satisfactory, we decided to consider data-based methods. The final goal is the fusion of model-based and data-based methods, as proposed in [17][18]. At first, we would like to check the efficiency of data-based methods to better isolate faults from attacks for simple cases.

The fault diagnosis and attack isolation problem can be formulated as a classical classification problem. The faults and attacks are considered as different classes in the classification problem. Thanks to simulation, we created a large set of data for each scenario and tried to apply common methods of the machine learning literature, focused on supervised learning with a classification goal, each class corresponding to a fault or an attack. These methods are basically Discriminant Analysis [START_REF] Friedman | Regularized discriminant analysis[END_REF], K-Nearest Neighbors [START_REF] Cover | Nearest neighbor pattern classification[END_REF], decision tree or bagged and boosted decision trees [START_REF] Breiman | Classification and Regression Trees[END_REF]. The goal was then to find a predictive model based on both input and output data.

Features

The two-tank system has mainly 6 observable variables that have been taken as features and that correspond to the sensor data for the state variables (water level in each tank) and to the inputs. The features are then

h m 1 , h m 2 , Q m p , U m b , U m 0 and U m
p . Note that the classification methods have been tested without taking into account the additional pressure sensors proposed for deriving residuals r 4 and r 5 .

Training Set

Each scenario has been run 10 times and each run varies because of noises included in the simulation. Each of them has a duration of 250s with a sampling period of 1s, so that each run corresponds to 250 input data. Each input data has been tagged by one of the 17 simplified tags, defined in Table 3. These simplified tags have been used because some scenarios seemed to be exactly the same physically. Future work could investigate the search of differences between them.

As all the scenarios already include some "normal" data (for example, from t=0s up to t=40s, all the data will be Diagnosability group

F S Physical explanation G 1 : f 1 [0001] T Failure pump G 2 : f 2 , f 3 , f 7 , f 8 , f 11 , f 12 , a 2 , a 3 , a 4 , a 10 [1100] T Impact on T 1 G 3 : f 4 , a 1
[1000] T Leakage in T 1 G 4 : f 5 , a 5 , a 6 , a 7 , a 8 , a 12 , a 13 , a 14

[0100] T Impact on T 2 G 5 : f 6

[0011] T Sensor fault for P 1 flow G 6 : f 9 , a 11

[1001] T Problem on Q p G 7 : f 10

[0010] T PI controller fault G 8 : f 0 , a 9

[0000] T Reply attack or nominal mode [110000] T Problem on V b G 3 : f 4 , a 1

[100000] T Leakage in T 1 G 41 : f 5 , a 5 , a 13 , a 14

[010000] T Impact on T 2 G 42 : a 6 , a 7 , a 8 , a 12

[010001] T Attack on T 2 G 5 : f 6

[001100] T Sensor fault for P 1 flow G 6 : f 9 , a 11

[100100] T Problem on Q p G 7 : f 10

[001000] T PI controller fault G 8 : f 0 , a 9

[000000] T Replay attack or nominal mode Table 3: Simplified tags for classification methods tagged f 0 a 9 ), we simulate the remaining 25 scenarios (from f 1 to f 12 and a 1 to a 8 ) over the 27 initial ones. A total number of 62500 tagged input data have then been stored.

The training was done with Matlab Classification Learner app, that performs supervised machine learning given a known set of input data and known responses to the data (here the tags).

Validation Scheme and Performances

The chosen validation scheme is the cross-validation method with 5 folds to partition the data set. The app partitions the data into 5 disjoint sets. For each fold, it trains a model using the out-of-fold observations, then assesses model performance using in-fold data.

Table 4 gives the performance results for the main methods included in the Matlab Classification Learner app.

The data-based methods have excellent results for classifying the different scenarios. Fine-kNN has an accuracy of 97.9% when only one neighbor is considered, while the accuracy is 97.4% when 5 neighbors are considered. The Bagged Trees method has an accuracy of 97.5%.

We show on Figure 3 and 4 the confusion matrices for the Fine-kNN and for the Bagged Trees methods. The confusion matrix gives the percentage of true class among the predicted class. As it is possible to see, the worst percentage of true detection is 78% for the Fine-kNN method and 82% for the Bagged Trees method.

From Figure 3 and 4, we can deduce that all the defined tagged cases can be detected by simple data-based methods. Faults and attacks can be isolated, except a 9 , a 1 and a 5 . These three attacks will be studied in future works. However, the isolation power of data-based methods is much better than model-based methods. For example, with modelbased methods, G 6 includes f 9 and a 11 while data-based methods succeed in isolate both cases. The case f 5 a 5 remains the most difficult case to isolate from others.

New Data Prediction

We exported the models deduced from the Fine-kNN method and from the Bagged-Tree method to the workspace of Matlab and used these trained models to make predictions using new data.

A new scenario involving a failure of the pump has been designed, with different flows. From 0 to 40s, the system has no failure and no attack. From 40s to 120s, the pump is faulty, then after 120s the system again in nominal mode. New data have been tested with the predicted models. Figure 5 shows the prediction results for the Fine-kNN method: each new data has been well classified, except during the transition from normal mode to faulty mode where one data is tagged f 6 instead of f 1 . The same problem occurs when using the model trained with Bagged Trees.

Another scenario has been tested, involving an attack corresponding of scenario a 13 . Figure 6 results for the Fine-kNN method: the attack is immediately detected and tagged a 13,14 with a Fine-KNN method, while the model-based methods kept the ambiguity with fault f 5 .

The same results are found with the Bagged Trees method.

Conclusion

This article presents the use of model-based diagnosis and data-based techniques for detection and isolation of cyberattacks and faults on a two-tank system. A set of 12 fault scenarios and 14 attack scenarios was described, 9 of which are new compared to previous works. A first model-based technique was successfully implemented and tested on the benchmark. All faults and attacks are detectable, except replay attack. By introducing 2 new residuals implemented on the benchmark the isolability has been increased so that some attacks can be isolated from faults. However, most of the time the different scenarios have the same signature. In this article the authors propose a first step towards fault and attack detection and isolation. More efficient methods have to be tested to take the fault and attacks dynamic into account. Model-based approaches in the field of switched systems have been mentioned. Another idea could be to take advantage of structural analysis to analyze internal dynamics of the model and residuals. Data-based diagnosis was performed using supervised machine learning methods. These methods obtained very good results in isolating faults and attacks and could be used in the future either alone, either using a fusion solution. Some aspects of the work however could be consolidated notably the characterization of the faults and attacks. Unsupervised learning methods were also tested. The goal is to find hidden pattern or groupings in data. This approach did not give satisfying results at the moment and should be investigated deeply. Methods for anomaly detection could be beneficial to consider cyberattacks not known in advance. Another research direction is data pre-processing to extract relevant features to fed machine learning methods. Here also dynamical aspects could be better considered.

Figure 2 :

 2 Figure 2: Residual values for scenario f 9

Figure 3 :Figure 4 :

 34 Figure 3: Confusion matrix for the Fine-kNN method

Figure 5 :

 5 Figure 5: Predicted classes on a scenario of type f 1

Figure 6 :

 6 Figure 6: Predicted classes on a scenario of type a 13

Table 1 :

 1 Diagnosability groups with 4 residuals Failure pump G 21 : f 2 , f 12 , a 2 , a 3 , a 4 , a 10 [110010] T Level sensor mh 1 involved G 22 : f 3 [110001] T Level sensor mh 2 involved G 23 : f 7 , f 8 , f 11

	Diagnosability group	F S	Physical explanation
	G 1 : f 1	[000100] T

Table 2 :

 2 Diagnosability groups with 6 residuals

	scenario tag	scenario	tag
	f 0 , a 9	f 0 a 9	f 9	f 9
	f 1	f 1	f 10	f 10
	f 2 , f 12	f 2,12 a 2 , a 3 , a 4 a 2,3,4
	f 3	f 3	a 6 , a 7 , a 8 a 6,7,8
	f 4 , a 1	f 4 a 1	a 10	a 10
	f 5 , a 5	f 5 a 5	a 11	a 11
	f 6	f 6	a 12	a 12
	f 7 , f 11	f 7,11 a 13 , a 14	a 13,14
	f 8	f 8

Table 4 :

 4 Performances Results for the main classification included in Matlab

	shows the prediction

Acknowledgments

We thank Julie Fouan, Maeva Kleinberg and Florian Valette, 3 th year students from INP-ENSEEIHT, Toulouse, who implemented the matlab simulator and the benchmark improvements presented in this paper.


Marta Marulli

Aurélie Edwards

Vuk Milišić

Nicolas Vauchelet

On the role of the epithelium in a model of sodium exchange in renal tubules

Keywords: Counter-current, transport equation, ionic exchange, stationary system, eigenproblem, long-time asymptotics

In this study we present a mathematical model describing the transport of sodium in a fluid circulating in a counter-current tubular architecture, which constitutes a simplified model of Henle's loop in a kidney nephron. The model explicitly takes into account the epithelial layer at the interface between the tubular lumen and the surrounding interstitium. In a specific range of parameters, we show that explicitly accounting for transport across the apical and basolateral membranes of epithelial cells, instead of assuming a single barrier, affects the axial concentration gradient, an essential determinant of the urinary concentrating capacity. We present the solution related to the stationary system, and we perform numerical simulations to understand the physiological behaviour of the system. We prove that when time grows large, our dynamic model converges towards the stationary system at an exponential rate. In order to prove rigorously this global asymptotic stability result, we study eigenproblems of an auxiliary linear operator and its dual.

Introduction

One of the main functions of the kidneys is to filter metabolic wastes and toxins from plasma and excrete them in urine. The kidneys also play a key role in regulating the balance of water and electrolytes, long-term blood pressure, as well as acid-base equilibrium. The structural and functional units of the kidney are called nephrons, which number about 1 million in each human kidney [START_REF] Bertram | Human nephron number: implications for health and disease[END_REF].

Blood is first filtered by glomerular capillaries and then the composition of the filtrate varies as it flows along different segments of the nephron : the proximal tubule, Henle's loop (which is formed by a descending limb and an ascending limb), the distal tubule, and the collecting duct. The reabsorption Figure 1: Simplified model of loop of Henle. q 1 /q 2 and u 1 /u 2 denote solute concentration in the epithelial layer and lumen of the descending/ascending limb, respectively. of water and solutes from the tubules into the surrounding interstitium (or secretion in the opposite direction) allows the kidneys to match precisely urinary excretion to the dietary intake [START_REF] Palmer | Integrated control of Na transport along the nephron[END_REF].

In the last decade, several groups have developed sophisticated models of water and electrolyte transport in the kidney. These models can be broadly divided into 2 categories: (a) detailed cell-based models that incorporate cellspecific transporters and predict the function of small populations of nephrons at steady-state ( [START_REF] Nieves-Gonzalez | Transport efficiency and workload distribution in a mathematical model of the thick ascending limb[END_REF]; [START_REF] Weinstein | A mathematical model of the rat nephron: glucose transport[END_REF]; [START_REF] Layton | A computational model for simulating solute transport and oxygen consumption along the nephrons[END_REF]; [START_REF] Weinstein | A mathematical model of the rat kidney: K + -induced natriuresis[END_REF]; [START_REF] Edwards | A model of uric acid transport in the rat proximal tubule[END_REF]), and (b) macroscale models that describe the integrated function of nephrons and renal blood vessels but without accounting for cell-specific transport mechanisms ( [START_REF] Tournus | Modèles d'échanges ioniques dans le rein: théorie, analyse asymptotique et applications numériques[END_REF], [START_REF] Tournus | Analysis of a simplified model of the urine concentration mechanism[END_REF]; [START_REF] Moss | Hormonal regulation of salt and water excretion: a mathematical model of whole kidney function and pressure natriuresis[END_REF]; [START_REF] Fry | Impact of nitric-oxide-mediated vasodilation and oxidative stress on renal medullary oxygenation: a modeling study[END_REF]; [START_REF] Clemmer | Mechanisms of blood pressure salt sensitivity: new insights from mathematical modeling[END_REF]; [START_REF] Hallow | A Quantitative Systems Physiology Model of Renal Function and Blood Pressure Regulation: Application in Salt-Sensitive Hypertension[END_REF]). These latter models do not consider explicitly the epithelial layer separating the tubule lumen from the surrounding interstitium, and represent the barrier as a single membrane. We developed the model presented below to assess the impact of this set of assumptions. Specifically, in this study we present a simplified mathematical model of solute transport in Henle's loop. This model accounts for ion transport between the lumen and the epithelial cells, and between the cells and the interstitium. The aim of this work is to evaluate the impact of explicitly considering the epithelium on predicted solute concentration gradients in the loop of Henle.

In our simplified approach, the loop of Henle is represented as two tubules in a counter-current arrangement, the descending and ascending limb are considered to be rigid cylinders of length L lined by a layer of epithelial cells. Water and solute reabsorption from the luminal fluid into the interstitium proceeds in two steps : water and solutes cross first the apical membrane at the lumencytosol interface and then the basolateral membrane at the cytosol-interstitium interface, [START_REF] Weinstein | A mathematical model of rat ascending Henle limb. II. Epithelial function[END_REF]. A schematic representation of the model is given in Figure 1.

The energy that drives tubular transport is provided by Na + /K + -ATPase, an enzyme that couples the hydrolysis of ATP to the pumping of sodium (Na + ) ions out of the cell and potassium (K + ) ions into the cell, across the basolateral membrane. The electrochemical potential gradients resulting from this active transport mechanism in turn drive the passive transport of ions across other transporters, via diffusion or coupled transport. We refer to diffusion as the biological process in which a substance tends to move from an area of high concentration to an area of low concentration [START_REF] Shultz | Basic Principles Of Membrane Transport[END_REF][START_REF] Thomas | Kidney modeling: status and perspectives[END_REF]. As described in [START_REF] Keener | Systems physiology[END_REF], in the absence of electrical forces, the diffusive solute flux from compartment 1 to compartment 2 (expressed in [mol.m -1 .s -1 ]) is given by:

J diffusion = P (u 2 -u 1 ),
where P [m.s -1 ] is the permeability of the membrane to the considered solute, the perimeter of the membrane, and u 1 and u 2 are the respective concentrations of the solute in compartments 1 and 2.

We assume that the volumetric flow rate in the luminal fluid (denoted by α > 0) remains constant, i.e. there is no transepithelial water transport. The descending limb is in fact permeable to water, but we make this simplifying assumption in order to facilitate the mathematical analysis. Given the countercurrent tubular architecture, this flow rate has a negative value in the ascending limb.

The present model focuses on tubular Na + transport. The concentration of Na + ([mol.m -3 ]) is denoted by u 1 and u 2 , respectively, in the lumen of the descending and ascending limb, by q 1 and q 2 in the epithelial cells of the descending and ascending limbs, respectively, and by u 0 in the interstitium. The permeability to Na + of the membrane separating the lumen and the epithelial cell of the descending and ascending limb is denoted by P 1 and P 2 , respectively. P 1,e denotes the permeability to Na + of the membrane separating the epithelial cell of the descending limb and the interstitium; the Na + permeability at the interface between the epithelial cell of the ascending limb and the interstitium is taken to be negligible. The re-absorption or secretion of ions generates electrical potential differences across membranes. In the present model, the impact of transmembrane potentials on Na + transport is not taken into account.

The concentrations depend on the time t and the spatial position x ∈ [0, L]. The dynamics of Na + concentration is given by the following model on (0, +∞)× (0, L)

a 1 ∂u 1 ∂t + α ∂u 1 ∂x = J 1 , a 2 ∂u 2 ∂t -α ∂u 2 ∂x = J 2 , (1)
a 3 ∂q 1 ∂t = J 3 , a 4 ∂q 2 ∂t = J 4 , a 0 ∂u 0 ∂t = J 0 . (2)
The parameters a i , for i = 0, 1, 2, 3, 4, denote positive constants defined as:

a 1 = πr 2 1 , a 2 = πr 2 2 , a 3 = π(r 2 1,e -r 2 1 ), a 4 = π(r 2 2,e -r 2 2 ), a 0 = π r 2 1,e + r 2 2,e 2 .
In these equations, r i , i = 1, 2, denotes the inner radius of tubule i, whereas r i,e denotes the outer radius of tubule i, which includes the epithelial layer. The fluxes J i describe the ionic exchanges between the different domains. They are modeled in the following way:

Lumen. In the lumen, we consider the diffusion of Na + towards the epithelium. Then,

J 1 = 2πr 1 P 1 (q 1 -u 1 ), J 2 = 2πr 2 P 2 (q 2 -u 2 ).
Epithelium. We take into account the diffusion of Na + from the descending limb epithelium towards both the lumen and the interstitium, J 3 = 2πr 1 P 1 (u 1 -q 1 ) + 2πr 1,e P 1,e (u 0 -q 1 ).

In the ascending limb (tubule 2), we also consider the active reabsorption that is mediated by Na + /K + -ATPase, which pumps 3 Na + ions out of the cell in exchange for 2 K + ions.

The net flux into the ascending limb epithelium is given by the sum of the diffusive flux from the lumen and the export across the pump, which is described using Michaelis-Menten kinetics, [START_REF] Keener | Systems physiology[END_REF]:

J 4 = 2πr 2 P 2 (u 2 -q 2 ) -2πr 2,e G(q 2 ),
where

G(q 2 ) = V m q 2 K M,2 + q 2 3 .
The exponent of G is related to the number of exchanged sodium ions. The affinity of the pump K M,2 , and its maximum velocity V m , are given real numbers. We notice that when q 2 → +∞, then G(q 2 ) → V m which is in accordance with the biological observation that the pump can be saturated.

Interstitium.

J 0 = 2πr 1 P 1,e (q 1 -u 0 ) + 2πr 2,e G(q 2 ).

We model the dynamics of a solute (here sodium) by the evolution of its concentration in each tubule. The transport of solute and its exchange are then modelled by a hyperbolic PDE system at constant speed with a non-linear transport term and with specific boundary conditions. The dynamics of ionic concentrations is given by the following model:

               a 1 ∂ t u 1 (t, x) + α∂ x u 1 (t, x) = J 1 (t, x) a 2 ∂ t u 2 (t, x) -α∂ x u 2 (t, x) = J 2 (t, x) a 3 ∂ t q 1 (t, x) = J 3 (t, x) a 4 ∂ t q 2 (t, x) = J 4 (t, x) a 0 ∂ t u 0 (t, x) = J 0 (t, x). (3)
We set the boundary conditions :

u 1 (t, 0) = u b (t), u 2 (t, L) = u 1 (t, L), t > 0, (4)
where u b is a given function in L ∞ (R + )∩L 1 loc (R + ), which is such that lim t→∞ u b (t) = ūb for some positive constant ūb > 0.

Finally, the system is complemented with initial conditions

u 1 (0, x) = u 0 1 (x), u 2 (0, x) = u 0 2 (x), u 0 (0, x) = u 0 0 (x), q 1 (0, x) = q 0 1 (x), q 2 (0, x) = q 0 2 (x).
To simplify notations in (3), we set K 1 := 2πr 1,e P 1,e , k 1 := 2πr 1 P 1 , and k 2 := 2πr 2 P 2 . For the diffusive fluxes J 4 and J 0 , we include the constant 2πr 2,e in the parameter V m and replace the parameter V m with V m,2 := 2πr 2,e V m , such that

G(q 2 ) = V m,2 q 2 K M,2 + q 2 3 . (5)
Moreover, the orders of magnitude of k 1 , k 2 are the same even if their values are not definitely equal, we may assume to further simplify the analysis that k 1 = k 2 = k. We will refer to this as the dynamic system and then (3) reads :

a 1 ∂ t u 1 + α∂ x u 1 = k(q 1 -u 1 ) (6a)
a 2 ∂ t u 2 -α∂ x u 2 = k(q 2 -u 2 ) (6b)
a 3 ∂ t q 1 = k(u 1 -q 1 ) + K 1 (u 0 -q 1 ) (6c)
a 4 ∂ t q 2 = k(u 2 -q 2 ) -G(q 2 ) (6d) a 0 ∂ t u 0 = K 1 (q 1 -u 0 ) + G(q 2 ). (6e)
The existence and uniqueness of vector solution u = (u 1 , u 2 , q 1 , q 2 , u 0 ) to this system are investigated in [15]. Several previous works have neglected the epithelium region (see e.g. [START_REF] Tournus | Analysis of a simplified model of the urine concentration mechanism[END_REF][START_REF] Tournus | Modèles d'échanges ioniques dans le rein: théorie, analyse asymptotique et applications numériques[END_REF]). The first goal of this work is to study the effects of this region in the mathematical model. The main indicator quantifying these effects is the parameter k which accounts for the permeability between the lumen and the epithelium. Then, we analyse the dependency between the concentrations and k. In the absence of physiological perturbations, the concentrations are very close to the steady state, thus it seems reasonable to consider solutions of (6) at equilibrium, which leads us to study the system :

                   +α∂ x ū1 = k(q 1 -ū1 ) -α∂ x ū2 = k(q 2 -ū2 ) 0 = k(ū 1 -q1 ) + K 1 (ū 0 -q1 ) 0 = k(q 2 -ū2 ) -G(q 2 ) 0 = K 1 (q 1 -ū0 ) + G(q 2 ) ū1 (L) = ū2 (L), ū1 (0) = ūb . (7)
Section 2 concerns the analysis of solutions to stationary system [START_REF] Fry | Impact of nitric-oxide-mediated vasodilation and oxidative stress on renal medullary oxygenation: a modeling study[END_REF]. In particular, we study their qualitative behaviour and their dependency with respect to the parameter k. Our mathematical observations are illustrated by some numerical computations.

The second aim of the paper is to study the asymptotic behaviour of the solutions of [START_REF] Feher | Quantitative Human Physiology[END_REF]. In Theorem 3.1, we show that they converge as t goes to +∞ to the steady state solutions solving [START_REF] Fry | Impact of nitric-oxide-mediated vasodilation and oxidative stress on renal medullary oxygenation: a modeling study[END_REF]. Section 3 is devoted to the statement and the proof of this convergence result. Finally, an Appendix provides some useful technical lemmas.

Stationary system

In this section, after proving basic existence and uniqueness results, we investigate how solutions of (7) depend upon the parameter k. We recall that it includes also the permeability parameter as k = k i with k i := 2πr i P i , i = 1, 2.

In order to study the qualitative behaviour of these solutions, we then perform some numerical simulations.

Stationary solution

We first show existence and uniqueness of solutions to the stationary system: Lemma 2.1. Let ūb > 0. Let G be a C 2 function, uniformly Lipschitz, such that G is uniformly bounded and G(0) = 0 (e.g. the function defined in (5)). Then, there exists an unique vector solution to the stationary problem [START_REF] Fry | Impact of nitric-oxide-mediated vasodilation and oxidative stress on renal medullary oxygenation: a modeling study[END_REF].

Moreover, if we assume that G > 0 on R + , then we have the following relation q2 < ū < q1 < ū0 .

Proof. Summing up all the equations of system ( 7), we deduce that α(∂ x ū1 -∂ x ū2 ) = 0. From the boundary condition ū1 (L) = ū2 (L), we obtain ū1 = ū2 = ū. Therefore, we may simplify system [START_REF] Fry | Impact of nitric-oxide-mediated vasodilation and oxidative stress on renal medullary oxygenation: a modeling study[END_REF] in

               α∂ x ū = k(ū -q2 ) 2ū = q1 + q2 0 = k(ū -q1 ) + K 1 (ū 0 -q1 ) 0 = k(ū -q2 ) -G(q 2 ) 0 = K 1 (q 1 -ū0 ) + G(q 2 ). (8)
By the fourth equation of ( 8), ū = q2 + G(q2) k , inserted into the first equation, it gives ∂ x ū = G(q2) α . We obtain a differential equation satisfied by q2 ,

∂ x q2 = G(q 2 ) α + α k G (q 2 ) , (9)
with α, k positive constants and provided with the initial condition q2 (0) that satisfies

q2 (0) + G(q 2 (0)) k = ūb . (10)
We first remark that q2 (0) → q2 (0) + G(q2(0)) k is a C 2 increasing function which takes the value 0 at 0 and goes to +∞ at +∞. Thus, for any ūb > 0 there exists a unique q2 (0) > 0 solving [START_REF] Keener | Systems physiology[END_REF].

By assumption, G and G are uniformly bounded, thus we check easily that the right-hand side of ( 9) is uniformly Lipschitz. Therefore, the Cauchy problem ( 9)-( 10) admits a unique solution, which is positive (by uniqueness since 0 is a solution).

Then, other quantities are computed thanks to the relations:

ū = q2 + G(q 2 ) k , q1 = q2 + 2G(q 2 ) k , ū0 = 1 K 1 + 2 k G(q 2 ) + q2 . (11)
Moreover, by the fourth and fifth equations of system (8) and since G(q 2 ) > 0, we immediately deduce that q2 < ū and q1 < ū0 . Using the second equation of ( 8), we obtain the claim.

Parameters

Description Values L Length of tubules 2 • 10 -3 [m] α
Water flow in the tubules

10 -13 [m 3 /s] r i Radius of tubule i = 1, 2 10 -5 [m] r i,e
Radius of epithelium layer i = 1, 2 1.5

• 10 -5 [m] K 1 2πr 1,e P 1,e ∼ 2π • 10 -11 [m 2 /s] k = k i 2πr i P i , i = 1, 2 changeable [m 2 /s] V m,2
Rate of active transport

∼ 2πr 2,e 10 -5 [mol.m -1 .s -1 ] K M,2
Pump affinity for sodium (N a + )

3, 5 [mol/m 3 ] ūb Initial concentration in tubule 1 140 [mol/m 3 ]
Table 1: Frequently used parameters

Numerical simulations of stationary solutions

We approximate numerically solutions of [START_REF] Hallow | A Quantitative Systems Physiology Model of Renal Function and Blood Pressure Regulation: Application in Salt-Sensitive Hypertension[END_REF]. Numerical values of the parameters (cf Table 1) are extracted from Table 2 in [START_REF] Hervy | Inner medullary lactate production and urineconcentrating mechanism: a flat medullary model[END_REF] and Table 1 in [START_REF] Layton | A region-based mathematical model of the urine concentrating mechanism in the rat outer medulla. I. Formulation and base-case results[END_REF].

Taking into account these quantities allow us to have the numerical ranges of the constants and the solution results in a biologically realistic framework. Following the proof of Lemma 2.1, we first solve [START_REF] Keener | Systems physiology[END_REF] thanks to a Newton method. Then, we solve ( 9) with a fourth order Runge-Kutta method. Finally, we deduce other concentrations u, q 1 , u 0 using [START_REF] Layton | Modeling transport and flow regulatory mechanisms of the kidney[END_REF].

Results from Figures 2a and2c show that in all compartments, concentrations increase as a function of depth (x-axis). Physiologically, this means that the fluid is more concentrated towards the hairpin turn (x = L) than near x = 0, because of active transport in the ascending limb. It can also be seen that Na + concentration is higher in the central layer of interstitium and lower in the ascending limb epithelium owing to active Na + transport from the latter to the central compartment, described by the non-linear term G(q 2 ). Furthermore, Figure 2b and Figure 2d highlight that increasing the permeability value homogenizes the concentrations in the tubules and in the epithelium region. Taking a very large permeability value is equivalent to fusing the epithelial layer with the adjacent lumen, such that luminal and epithelial concentrations become equal. It is proved rigorously in [15] that this occurs in the dynamic system [START_REF] Feher | Quantitative Human Physiology[END_REF]. This is derived and explained formally in Appendix A.

Figures 3 and4 depict the impact of permeability P 1 = P 2 = P on concentration profiles for various pump rates V m,2 . Axial profiles of luminal concentrations are shown in Figures 3a and4a, considering different values of the permeability between the lumen and the epithelium. The fractional increase in concentration (FIC) is shown in Figures 3b and4b : for each permeability value (plotted on the horizontal axis), we compute the following ratio (shown on the vertical axis):

FIC(ū) := 100 ū(L) -ū(0) ū(0) , (12)
where ū(L) is the concentration in the tubular lumen 1, 2 at x = L and ū(0) the concentration at x = 0. This illustrates the impact of permeability on the axial concentration gradient. We observe that this ratio depends also strongly on the value of V m,2 . The permeability range (numerically P ∈ [10 -8 , 10 -5 ], equispaced 50 values between these) encompasses the physiological value which should be around 10 -7 m/s. As shown in Figures 3b and4b, the FIC increases significantly with P until it reaches a plateau : indeed, as diffusion becomes more rapid than active transport (that is, pumping by Na + / K + -ATPase), the permeability ceases to be rate-limiting. As shown by comparing Figures 3b and4b, the FIC is strongly determined by the pump rate V m,2 : if P ∈ [10 -8 , 10 -6 ], Na + concentration along the lumen increases by less than 12% if V m,2 = 2πr 2,e 10 -5 , and may reach 120% if V m,2 = 2πr 2,e 10 -4 . This raise is expected since concentration differences are generated by active transport ; the higher the rate of active transport, the more significant these differences. Conversely, in the absence of pumping, concentrations would equilibrate everywhere. In regards to the axial gradient, the interesting numerical results are in Figure (4a) and (4b). We observe that the axial gradient increases with increasing permeability when the latter is varied within the chosen range. Therefore this indicates that taking into account the epithelial layer in the model has a significant influence on the axial concentration gradient.

Moreover, numerical results also confirm that : ū1 = ū2 = u < q1 < ū0 as reported in Lemma 2.1. We recall that we assume a constant water flow α which allows us to deduce ū1 = ū2 . As noted above, the descending limb is in fact very permeable to water and α should decrease significantly in this tubule, such that ū1 differs from ū2 , except at the hairpin turn at x = L. On the other hand, the last equation of system [START_REF] Hallow | A Quantitative Systems Physiology Model of Renal Function and Blood Pressure Regulation: Application in Salt-Sensitive Hypertension[END_REF] implies that q1 < ū0 , meaning that the concentration of Na + is lower in the epithelial cell than in the interstitium, as observed in vivo, [START_REF] Atherton | Time course of changes in renal tissue and urinary composition after cessation of constant infusion of lysine vasopressin in the conscious, hydrated rat[END_REF].

With the expression of G in (5), equation [START_REF] Keener | Systems physiology[END_REF] reads

q2 (0) + V m,2 k q2 (0) K M,2 + q2 (0) 3 = ūb . (13)
In order to better understand the behaviour of the axial concentration gradient shown in Figures (3b), (4b), we compute the derivative of ( 13) with respect to the parameter V m and with respect to k respectively:

∂ q2 (0) ∂V m + 1 k G (q 2 (0)) ∂ q2 (0) ∂V m + 1 k q2 (0) K M,2 + q2 (0) 3 = 0, ∂ q2 (0) ∂k + 1 k G (q 2 (0)) ∂ q2 (0) ∂k - 1 k 2 G(q 2 (0)) = 0. Then, we get ∂ q2 (0) ∂V m = -1 k 1 + 1 k G (q 2 (0)) q2 (0) k M + q2 (0) 3 ≤ 0, ∂ q2 (0) ∂k = 1 k 2 G(q 2 (0)) 1 + 1 k G (q 2 (0)) ≥ 0,
because G is a monotone non-decreasing function and q 2 (0) is positive.

We observe from numerical results (see Figures 2a,2c, 3a, and 4a) that the gradient of u is almost constant. Thus, we may make the approximation

∂ x ū ∼ ∂ x ū(0) = G(q 2 (0)) α . (14)
Its derivatives with respect to V m,2 and k are both non negative :

∂ ∂k [∂ x ū] ∼ G (q 2 (0)) α G(q2(0)) k 2 1 + 1 k G (q 2 (0)) ≥ 0, ∂ ∂V m [∂ x ū] ∼ 1 α q2 (0) K M + q2 (0) 3 1 1 + 1 k G (q 2 (0)) ≥ 0.
It means that the axial concentration gradient is an increasing function both with respect to the rate of active transport V m,2 and to the permeability k. Indeed in Fig. 5, we perform numerical simulations varying both P and V m,2 and observe that the ratio (12) increases monotonically with respect to both parameters.

Limiting cases: k → ∞, and k → 0

Numerical results show that above a certain high value of permeability, the epithelial concentration in tubule 2 seems to reach a plateau (see Figures 3b and4b). There are two different regimes : one for large values of permeabilities, one for small values of permeabilities, and a fast transition between them.

In the large permeabilities asymptotic, we may approximate system (8) by the limiting model k = +∞. In this case, [START_REF] Layton | Modeling transport and flow regulatory mechanisms of the kidney[END_REF]

reduces to ū = q2 , q1 = q2 , ū0 = G(q 2 ) K 1 + q2 .
for all x ∈ (0, L). This is understandable from a formal point of view, also taking into account computations in Appendix (A) for the stationary system [START_REF] Fry | Impact of nitric-oxide-mediated vasodilation and oxidative stress on renal medullary oxygenation: a modeling study[END_REF]. In this case, the gradient concentration is directly proportional to V m,2 :

∂ x q2 = G(q 2 ) α + α k G (q 2 ) -→ k→+∞ G(q 2 ) α = ∂ x ū.
From ( 9), the Cauchy problem reduces to

∂ x q2 (x) = G(q 2 ) α , q2 (0) = ūb .
Additionally, it is clear that the higher pump value, the more the FIC will increase, as observed in Figure 4b.

On the other hand for small values of permeability, we obtain formally

∂ x q2 -→ k→0 0, ∂ x ū = G(q 2 ) α .
Therefore, in a neighbourhood of the value P ∼ 10 -8 , the concentration gradient tends to be constant and for this reason we notice a plateau.

Long time behaviour

This section is devoted to the main mathematical result of this paper concerning the long time asymptotics of solutions to (6) towards solutions to the stationary system [START_REF] Fry | Impact of nitric-oxide-mediated vasodilation and oxidative stress on renal medullary oxygenation: a modeling study[END_REF] as time goes to +∞. We first state the main result and the assumptions needed. Then, we introduce eigenelements of an auxiliary linear system and its dual problem. Using these auxiliary functions, we are able to show the convergence when the time variable goes to +∞. A similar approach was considered in [START_REF] Tournus | Analysis of a simplified model of the urine concentration mechanism[END_REF] following ideas from [START_REF] Perthame | Exponential decay for the fragmentation or cell-division equation[END_REF].

Statement of the main result

Before stating the main result, we provide assumptions on the initial and boundary data.

Assumption 3.1. We assume that the initial solute concentrations are nonnegative and uniformly bounded in L ∞ (0, L) and in the total variation :

0 ≤ u 0 1 , u 0 2 , q 0 1 , q 0 2 , u 0 0 ∈ BV (0, L) ∩ L ∞ (0, L). ( 15
)
Assumption 3.2. The boundary condition of system (6) is such that

0 ≤ u b ∈ L ∞ (R + ) ∩ L 1 loc (R + ), lim t→+∞ |u b -ūb | = 0, (16)
for some constant ūb > 0.

BV is the space of functions with bounded variation, we notice that such functions have a trace on the boundary (see e.g. [START_REF] Evans | Measure Theory and Fine Properties of Functions[END_REF]); hence the boundary condition u 2 (t, L) = u 1 (t, L) is well-defined. Assumption 3.3. Regularity and boundedness of G. We assume that the nonlinear function modelling active transport in the ascending limb (tube 2) is a bounded and Lipschitz-continuous function on R + :

∀x ∈ R + , 0 ≤ G(q 2 ) ≤ G ∞ , 0 ≤ G (q 2 ) ≤ G ∞ . (17)
We notice that G defined by ( 5) satisfies straightforwardly [START_REF] Nieves-Gonzalez | Transport efficiency and workload distribution in a mathematical model of the thick ascending limb[END_REF]. We now state the main result.

Theorem 3.1 (Long time behaviour). Under Assumptions 3.1, 3.2 and 3.3, the solution to the dynamical problem (6) denoted by u(t, x) = (u 1 , u 2 , q 1 , q 2 , u 0 ) converges as time t goes to +∞ towards ū(x), the unique solution to the stationary problem [START_REF] Fry | Impact of nitric-oxide-mediated vasodilation and oxidative stress on renal medullary oxygenation: a modeling study[END_REF], in the following sense

lim t→+∞ u(t) -ū L 1 (Φ) = 0,
with the space

L 1 (Φ) = u : [0, L] → R 5 ; u L 1 (Φ) := L 0 |u(x)| • Φ(x) dx < ∞ ,
where Φ = (ϕ 1 , ϕ 2 , φ 1 , φ 2 , ϕ 0 ) is defined in Proposition 3.1 below. Moreover, if we assume that there exist µ 0 > 0 and C 0 such that |u b (t)-ū b | ≤ C 0 e -µ0t for all t > 0, then there exist µ > 0 and C > 0 such that we have the convergence with an exponential rate

u(t) -ū L 1 (Φ) ≤ Ce -µt . ( 18
)
The scalar product used in the latter claim means :

L 0 |u(x)| • Φ(x) dx = L 0 |u 1 |ϕ 1 (x) + |u 2 |ϕ 2 (x) + |q 1 |φ 1 (x) + |q 2 |φ 2 (x) + |u 0 |ϕ 0 (x) dx.
The definition of the left eigenvector Φ and its role are given hereafter.

The eigen-problem

In order to study the long time asymptotics of the time dependent system (6), we consider the eigen-problem associated with a specific linear system [START_REF] Perthame | Transport equations in biology[END_REF][START_REF] Tournus | Analysis of a simplified model of the urine concentration mechanism[END_REF]. This system is, in some sort, a linearized version of the stationary system [START_REF] Fry | Impact of nitric-oxide-mediated vasodilation and oxidative stress on renal medullary oxygenation: a modeling study[END_REF] where the derivative of the non-linearity is replaced by a constant g. When these eigenelements (λ, U, Φ) exist, the asymptotic growth rate in time for a solution u of ( 6) is given by the first positive eigenvalue λ and the asymptotic shape is given by the corresponding eigenfunction U.

Let us introduce the eigenelements of an auxiliary stationary linear system

               ∂ x U 1 = λU 1 + k(Q 1 -U 1 ) -∂ x U 2 = λU 2 + k(Q 2 -U 2 ) 0 = λQ 1 + k(U 1 -Q 1 ) + K 1 (U 0 -Q 1 ) 0 = λQ 2 + k(U 2 -Q 2 ) -gQ 2 0 = λU 0 + K 1 (Q 1 -U 0 ) + gQ 2 , ( 19
)
where g is a positive constant which will be fixed later. This system is complemented with boundary and normalization conditions :

U 1 (0) = 0, U 1 (L) = U 2 (L), L 0 (U 1 + U 2 + Q 1 + Q 2 + U 0 ) dx = 1. ( 20
)
We also consider the related dual system :

               -∂ x ϕ 1 = λϕ 1 + k(φ 1 -ϕ 1 ) ∂ x ϕ 2 = λϕ 2 + k(φ 2 -ϕ 2 ) 0 = λφ 1 + k(ϕ 1 -φ 1 ) + K 1 (ϕ 0 -φ 1 ) 0 = λφ 2 + k(ϕ 2 -φ 2 ) + g(ϕ 0 -φ 2 ) 0 = λϕ 0 + K 1 (φ 1 -ϕ 0 ), (21)
with following conditions :

ϕ 1 (L) = ϕ 2 (L), ϕ 2 (0) = 0, L 0 (U 1 ϕ 1 + U 2 ϕ 2 + Q 1 φ 1 + Q 2 φ 2 + U 0 ϕ 0 ) dx = 1.
(22) For a given λ, the function U := (U 1 , U 2 , Q 1 , Q 2 , U 0 ) is the right eigenvector solving [START_REF] Perthame | Exponential decay for the fragmentation or cell-division equation[END_REF], while Φ := (ϕ 1 , ϕ 2 , φ 1 , φ 2 , ϕ 0 ) is the left one, associated with the adjoint operator. The following result shows the existence of a positive eigenvalue and some properties of eigenelements. We underline that in order to make the proof easier, we consider the case k = k 1 = k 2 but the same result could be extended to the more general case where k 1 = k 2 . Proposition 3.1. Let g > 0 be a constant. There exists a unique (λ, U, Φ) with λ ∈ (0, λ -) solution to the eigenproblem (19)-( 22), where

λ -= (2K 1 + k) -4K 2 1 + k 2 2
.

Moreover, we have U(x) > 0, Φ(x) > 0 on (0, L) and φ 2 < ϕ 0 .

In order to prove this result, we will divide the proof in two steps : Lemmas 3.1 and 3.2 respectively. Proposition 3.1 is a direct consequence of these two Lemmas. We start with the direct problem : Lemma 3.1 (The direct problem). There exists a unique λ > 0 such that the direct problem (19)-( 20) admits a unique positive solution U = (U 1 , U 2 , Q 1 , Q 2 , U 0 ) on (0, L), and 0 < λ < λ -.

Proof. Summing all equations in [START_REF] Perthame | Exponential decay for the fragmentation or cell-division equation[END_REF] we find that :

U 1 -U 2 = λ(U 1 + U 2 + Q 1 + Q 2 + U 0 ). ( 23
)
Integrating with respect to x and using condition [START_REF] Perthame | Transport equations in biology[END_REF], we obtain U 2 (0) = λ.

By the fourth equation in [START_REF] Perthame | Exponential decay for the fragmentation or cell-division equation[END_REF], we find directly:

Q 2 (x) = kU 2 (x) k + g -λ = U 2 (x) 1 + 1 k (g -λ) . ( 24
)
Putting this expression into the second equation in [START_REF] Perthame | Exponential decay for the fragmentation or cell-division equation[END_REF], we find

-U 2 = U 2 λ + λ -g 1 + 1 k (g -λ)
.

Solving the latter equation, we deduce that

U 2 (x) = U 2 (0)e -λx+ x 0 -λ+g 1+ 1 k (g-λ) dy = λe (-λ+η(λ))x ; with η(λ) := -λ + g 1 + 1 k (g -λ) . ( 25
)
Using the fifth equation of system [START_REF] Perthame | Exponential decay for the fragmentation or cell-division equation[END_REF] we recover

U 0 (x) = K 1 K 1 -λ Q 1 (x) + g K 1 -λ Q 2 (x).
We inject this into the third equation to obtain

Q 1 (x) k -λ - K 1 λ K 1 -λ = gK 1 K 1 -λ Q 2 (x) + kU 1 (x).
Thanks to [START_REF] Tournus | Modèles d'échanges ioniques dans le rein: théorie, analyse asymptotique et applications numériques[END_REF] we write also:

Q 1 (x) k -λ - K 1 λ K 1 -λ = K 1 g K 1 -λ 1 (1 + 1 k (g -λ)) U 2 (x) + kU 1 (x).
Taking into account the first equation of system [START_REF] Perthame | Exponential decay for the fragmentation or cell-division equation[END_REF], we obtain :

U 1 (x) = c λ U 1 (x) + k λ g 1 + 1 k (g -λ) U 2 (x), (26)
where we simplify notations by introducing :

k λ := k K1 K1-λ k -λ -K1λ K1-λ , c λ := λ + k(λ + K1λ K1-λ ) k -λ -K1λ K1-λ . ( 27
)
The denominator k -λ -K1λ K1-λ vanishes for

λ ± = (2K 1 + k) ± 4K 2 1 + k 2 2 .
Obviously lim λ→λ-k λ = +∞ and we also have that 0 < λ -< min(K 1 , k).

Now we solve directly the ODE [START_REF] Tournus | A model of calcium transport along the rat nephron[END_REF] with its initial condition, we get

U 1 (x) = λgk λ 1 + 1 k (g -λ)
e c λ x -e (η(λ)-λ)x c λ + λ -η(λ) .

We are looking for a λ > 0 such that boundary condition

U 1 (L) = U 2 (L) is satisfied, in other words U1(L) U2(L) = 1, namely F (λ) := gk λ 1 + 1 k (g -λ) e (c λ +λ-η(λ))L -1 c λ + λ -η(λ) = 1, (28)
where we recall that k λ , c λ are defined in [START_REF] Weinstein | A mathematical model of the rat nephron: glucose transport[END_REF] and η(λ) in [START_REF] Tournus | Analysis of a simplified model of the urine concentration mechanism[END_REF]. We remark immediately that for λ = 0 in ( 27), we have k 0 = 1, c 0 = 0. Then,

F (0) = 1 -exp - gL 1 + g k < 1.
We notice that for k λ , c λ > 0, F (λ) is a continuous increasing function with respect to λ since the product of increasing and positive functions is still increasing (see Appendix (B.1) for more details). Moreover lim λ→λ-F (λ) = +∞. Then it exists a unique λ ∈ (0, λ -) such that

F (λ) = 1. Moreover, for 0 < λ < λ -< min(k, K 1 ), the functions U 1 , U 2 , Q 1 , Q 2 , U 0 are positive on [0, L].
Lemma 3.2 (The dual problem). Let λ and U be as in Lemma (3.1). Then, there exists Φ := (ϕ 1 , ϕ 2 , φ 1 , φ 2 , ϕ 0 ), the unique solution of dual problem (21)-( 22) with ϕ 1 , ϕ 2 , φ 1 , φ 2 , ϕ 0 > 0. Moreover, we have φ 2 < ϕ 0 .

Proof. By the fifth equation of system [START_REF] Sands | The thick ascending limb and water channels: half-full or half-empty[END_REF] we have directly :

ϕ 0 = K 1 K 1 -λ φ 1 .
Replacing this expression in the third equation we obtain

(k -λ - K 1 λ K 1 -λ )φ 1 = kϕ 1 .
Then,

ϕ 0 (x) = k K1 K1-λ k -λ -K1λ K1-λ ϕ 1 (x) = k λ ϕ 1 (x),
where k λ is defined in [START_REF] Weinstein | A mathematical model of the rat nephron: glucose transport[END_REF]. Using the first equation of ( 21), we have

-ϕ 1 = ϕ 1 λ + k λ + K1λ K1-λ k -λ -K1λ K1-λ
.

Integrating, we obtain

ϕ 1 (x) = ϕ 1 (0)e -λx e -βx , β = β λ = λk(2K 1 -λ) λ 2 -2K 1 λ -λk + K 1 k . ( 29
)
We easily check that

β > 0 if 0 < λ < λ -< K 1 .
As shown in details in the Appendix B.2, for all x ∈ (0, L), (U

1 ϕ 1 ) - (U 2 ϕ 2 ) = 0. Integrating, we get U 1 (x)ϕ 1 (x) -U 2 (x)ϕ 2 (x) = U 1 (0)ϕ 1 (0) - U 2 (0)ϕ 2 (0) = 0, thanks to boundary conditions U 1 (0) = 0 and ϕ 2 (0) = 0.
(Notice also that taking x = L in this latter relation, and using the boundary condition U 1 (L) = U 2 (L) = 0, we recover ϕ 1 (L) = ϕ 2 (L).) Therefore, we get

ϕ 2 (x) = U 1 (x) U 2 (x) ϕ 1 (x), ∀x ∈ [0, L]. (30)
Using the fourth equation in [START_REF] Sands | The thick ascending limb and water channels: half-full or half-empty[END_REF] and thanks to (30), we obtain

0 = λφ 2 + k U 1 U 2 ϕ 1 -kφ 2 + gϕ 1 K 1 K 1 -λ • k k -λ -K1λ K1-λ
, which allows to compute φ 2 :

φ 2 (x) = k U1(x) U2(x) ϕ 1 (x) k -λ + g + g k -λ + g k λ ϕ 1 (x).
Each function depends on the first component of Φ, i.e. ϕ 1 (x), and to sum up, the following relation has been obtained:

                   ϕ 1 (x) = ϕ 1 (0)e -λx e -βx ϕ 2 (x) = U1(x) U2(x) ϕ 1 (x) φ 1 (x) = ϕ 1 (x) k k-λ- K 1 λ K 1 -λ φ 2 (x) = 1 k-λ+g k U1(x) U2(x) + gk λ ϕ 1 (x) ϕ 0 (x) = k λ ϕ 1 (x), (31)
where k λ , β are defined in [START_REF] Weinstein | A mathematical model of the rat nephron: glucose transport[END_REF] and [START_REF] Weinstein | A mathematical model of rat ascending Henle limb. II. Epithelial function[END_REF]. Hence the sign of Φ depends on the sign of ϕ 1 (0), the other quantities and constants being positive for λ ∈ (0, λ -) and g > 0 by assumption. Then, we use the normalization condition ( 22) and (31) in order to show the positivity of ϕ 1 (0). It implies that

ϕ 1 (0) L 0 e -λx e -βx 2U 1 (x) + Q 1 (x) k k -λ -K1λ K1-λ + + Q 2 (x) k + g + λ k U 1 (x) U 2 (x) + gk λ + k λ U 0 (x) dx = 1.
The integral on the left hand side is positive, thanks to properties of functions previously defined. Given that ϕ 1 (0) is constant and all other quantities positive, we can conclude that ϕ 1 (0) > 0.

We are left to prove that the quantity φ 2 -ϕ 0 is negative. Using (31), we rewrite :

φ 2 -ϕ 0 = kk λ ϕ 1 (x) k -λ + g 1 k λ U 1 U 2 + λ k -1 .
From the explicit expression of U 1 and U 2 , we have

φ 2 -ϕ 0 = k λ ϕ 1 (x) 1 + 1 k (g -λ)   x 0 g 1+ 1 k (g-λ) e -c λ (y-x) e (-λ+η(λ))y dy e (-λ+η(λ))x + λ k -1   = k λ ϕ 1 (x) 1 + 1 k (g -λ) g 1 + 1 k (g -λ) 1 -e -c λ x-λx+η(λ)x c λ + λ -η + λ k -1 ,
where we recall the notation η(λ) = -λ+g

1+ 1 k (g-λ)
. We set

H(x) := g 1 + 1 k (g -λ) 1 -e -c λ x-λx+η(λ)x c λ + λ -η(λ) . ( 32
)
We have

φ 2 -ϕ 0 < 0 ⇐⇒ H(x) + λ k -1 < 0.
We observe that H(0) = 0 and H(L) = 1 k λ thanks to [START_REF] Weinstein | A mathematical model of the rat kidney: K + -induced natriuresis[END_REF]. Moreover, H(L) < 1 -λ k for λ ∈ (0, λ -). Indeed, we have

λ < k - k k λ ,
which holds if and only if λ 2 -K1λ-kλ

K1

< 0 which is true on (0, λ -), since λ -< k by definition. Moreover, it is clear that H is an increasing function on [0, L] for λ ∈ (0, λ -). Then H(x) ≤ H(L) < 1 -λ k . This concludes the proof.

Proof of Theorem 3.1

Now we are ready to prove Theorem 3.1. We set d i (t, x) := |u i (t, x) -ūi (x)| i = 0, 1, 2 and δ j := |q j (t, x) -qj (x)|, j = 1, 2 with ūi , qi satisfying (7) and u i , q i solving (6). We subtract component-wise ( 6) to [START_REF] Fry | Impact of nitric-oxide-mediated vasodilation and oxidative stress on renal medullary oxygenation: a modeling study[END_REF]. Then we multiply each of the entries by sign(u i -ūi ) or sign(q j -qj ) respectively. We obtain the following inequalities :

               a 1 ∂ t d 1 + α∂ x d 1 ≤ k(δ 1 -d 1 ) a 2 ∂ t d 2 -α∂ x d 2 ≤ k(δ 2 -d 2 ) a 3 ∂ t δ 1 ≤ k(d 1 -δ 1 ) + K 1 (d 0 -δ 1 ) a 4 ∂ t δ 2 ≤ k(δ 2 -d 2 ) - Ĝ a 0 ∂ t d 0 ≤ K 1 (δ 1 -d 0 ) + Ĝ, (33)
with Ĝ := |G(q 2 ) -G(q 2 )|. We have used also the monotonicity of G (see [START_REF] Nieves-Gonzalez | Transport efficiency and workload distribution in a mathematical model of the thick ascending limb[END_REF]). We set

M (t) := L 0 (a 1 d 1 ϕ 1 + a 2 d 2 ϕ 2 + a 3 δ 1 φ 1 + a 4 δ 2 φ 2 + a 0 d 0 ϕ 0 ) dx.
Multiplying each equation of (33) by the corresponding dual function ϕ i , φ i , adding all equations and integrating with respect to x, we obtain :

d dt M (t) ≤ L 0 k(δ 1 -d 1 )ϕ 1 + k(δ 2 -d 2 )ϕ 2 + k(d 1 -δ 1 )φ 1 + K 1 (d 0 -δ 1 )φ 1 +k(δ 2 -d 2 )φ 2 -Ĝφ 2 +K 1 (δ 1 -d 0 )ϕ 0 + Ĝϕ 0 dx+α L 0 (∂ x d 2 ϕ 2 -∂ x d 1 ϕ 1 ) dx.
Integrating by parts the last integral and using the dual system [START_REF] Sands | The thick ascending limb and water channels: half-full or half-empty[END_REF], we can simplify the latter inequality into

d dt M (t) ≤ -λ L 0 (d 1 ϕ 1 + d 2 ϕ 2 + δ 1 φ 1 + δ 2 φ 2 + d 0 ϕ 0 ) dx +d 2 (L)ϕ 2 (L)-d 1 (L)ϕ 1 (L)-d 2 (0)ϕ 2 (0)+d 1 (0)ϕ 1 (0)+ L 0 (gδ 2 -Ĝ)(φ 2 -ϕ 0 ) dx.
Using the normalization conditions in [START_REF] Perthame | Transport equations in biology[END_REF] and in [START_REF] Shultz | Basic Principles Of Membrane Transport[END_REF], we obtain

d dt M (t) ≤ -λ max{a 1 , a 2 , a 3 , a 4 , a 0 } M (t)+d 1 (t, 0)ϕ 1 (0)+ L 0 (gδ 2 -Ĝ)(φ 2 -ϕ 0 ) dx.
To simplify the notation we set λ = -λ max{a1,a2,a3,a4,a0} . Since G is Lipschitzcontinuous and by assumption [START_REF] Nieves-Gonzalez | Transport efficiency and workload distribution in a mathematical model of the thick ascending limb[END_REF], Ĝ ≤ gδ 2 with g = G ∞ . With this choice of g, we apply Proposition 3.1 and deduce that the quantity

(φ 2 -ϕ 0 ) is negative. Then, d dt M (t) + λM (t) ≤ d 1 (t, 0)ϕ 1 (0).
Thanks to [START_REF] Moss | Hormonal regulation of salt and water excretion: a mathematical model of whole kidney function and pressure natriuresis[END_REF] and applying Gronwall's lemma, we conclude that M (t) ≤ M (0)e -λt + ϕ 1 (0) The first term of the right hand side is arbitrarily small at t goes to +∞. Hence, we have proved that for any ε > 0 there exists τ large enough such that for every t ≥ τ , M (t) ≤ M (0)e -λt + Cε.

Since M (t) = u(t) -ū L 1 (Φ) , it proves the convergence as stated in Theorem 3.1.

Finally, if we assume that there exist positive constants µ 0 and C 0 such that |u b (t) -ūb | ≤ C 0 e -µ0t , then from (34) we deduce

M (t) ≤ M (0)e -λt + C 0 ϕ 1 (0) e -µ0t -e -λt λ -µ 0 ≤ Ce -min{ λ,µ0}t .

Conclusion and outlook

In this study we present a model describing the transport of sodium in a simplified version of the loop of Henle in a kidney nephron. From a modelling point of view, it seems important to take into account the epithelium in the countercurrent tubular architecture since we observe that it may affect strongly the solute concentration profiles for a particular range of permeabilities.

The main limitation of the model is to not consider the re-absorption of water in descending limb. Indeed, in Section 2, we study the steady state solution and the assumption of a constant rate α and the boundary conditions lead to ū1 (x) = ū2 (x), i.e. the luminal concentrations of sodium are the same in both tubules for every x ∈ (0, L). Conversely, in vivo, the concentrations in lumen 1 and 2 are different due to the constitutive differences between the segments and presence of membrane channel proteins, for example the aquaporins. The thin descending limb of Henle's loop has low permeability to ions and urea, while being highly permeable to water. The thick ascending limb is impermeable to water, but it is permeable to ions. For this reason, a possible extension of the model shall assume that α is not constant but space-dependent. A first step could be, for instance, to take two different values of α for the first and second equation of the model ( 6), α 1 and α 2 . From the mathematical viewpoint, this choice slightly changes the structure of the hyperbolic system : for example, conservation of certain quantities should not be that easy to prove. Furthermore, this assumption about α has a relevant influence on other factors. As already pointed out, the relation between q1 and ū0 is biologically correct and consistent, this means that in vivo the concentration of Na + in the epithelial cell (intracellular) is lower than in interstitium. The intracellular concentrations (epithelium, q1 and q2 ) are usually of the order of 10mM whereas the extracellular ones (therefore in the lumen and in the interstitium) are of the order of 140mM, (see [START_REF] Feher | Quantitative Human Physiology[END_REF], page 692). There are also other types of source terms in the interstitium that could be added, accounting for blood vessels and/or collecting ducts. In this case, the last equation (6e) of the dynamic system should include a term that accounts for interstitium concentration storage or accumulation and for secretion-reabsorption of water and solutes, but the impact of adding such complex mechanisms in the model remains to be assessed.

In this study, we focused our attention on the axial concentration gradient and the FIC, previously defined in Section (2), which are significant factors in the urinary concentration mechanism, [START_REF] Layton | Modeling transport and flow regulatory mechanisms of the kidney[END_REF][START_REF] Layton | Mathematical Modeling in Renal Physiology[END_REF]. The axial gradient is an important determinant of urinary concentration capacity. When water intake is limited, mammals can conserve water in body fluids by excreting solutes in a reduced volume of water, that is, by producing a concentrated urine. The thick ascending limb plays an essential role in urine concentration and dilution, [START_REF] Sands | The thick ascending limb and water channels: half-full or half-empty[END_REF]: the active reabsorption of sodium without parallel reabsorption of water generates an interstitial concentration gradient in the outer medulla that in turn drives water reabsorption by the collecting ducts, thereby regulating the concentration of final urine. In summary, our model confirms that the active trans-epithelial transport of Na from the ascending limbs into the surrounding environment is able to generate an osmolality gradient. Our model indicates that explicitly accounting for the 2-step transport across the epithelium significantly impacts the axial concen-tration gradient within the physiological range of parameters values considered here. Thus, representing the epithelial layer as two membrane in series, as opposed to a single-barrier representation, may provide a more accurate understanding of the forces that contribute to the urinary concentrating mechanism.

B Technical results

B.1 Function F (λ)

In this subsection we prove the monotonicity of the function F (λ) which appears in the proof of Lemma 3.1. First let's recall it

F (λ) := gk λ 1 + 1 k (g -λ)
e (c λ +λ-η(λ))L -1

c λ + λ -η(λ) . ( 39
)
Lemma B.1. The function F defined by (39) is monotonically increasing on (0, λ -).

Proof. The product of positive increasing functions is increasing.

• λ → k λ = K1k λ 2 -2K1λ-kλ+kK1 is a positive and increasing function if λ ∈ (0, λ -). Indeed ∂k λ ∂λ = -2λkK1+2kK 2 1 +k 2 K1

(λ 2 -2K1λ-kλ+kK1) 2 is positive for 0 < λ < K 1 + k 2 and λ -< K 1 by definition.

• We set f 1 (λ) := g 1+ 1 k (g-λ) ; if λ < g +k the function f is positive since g > 0 by hypothesis and it is also increasing since

∂ ∂ λ f 1 (λ) = g k (1+ 1
k (g(y)-λ)) 2 > 0, and λ -≤ k 2 . • The function x → e x -1

x is increasing on R + and the function λ → c λ + λ -η(λ) is increasing on (0, λ 1 ). Indeed, we have straightforwardly

c λ + λ -η(λ) = 2λ + 2k + k 2 k -λ -K1λ K1-λ + k 2 k + g -λ .

B.2 Relation between direct and dual system

We recall the eigenelements problem written as below:

      ∂ x U 1 (x) -∂ x U 2 (x) 0 0 0       = λU(x) + AU(x); U(x) =       U 1 U 2 Q 1 Q 2 U 0       (40)       -∂ x ϕ 1 (x) ∂ x ϕ 2 (x) 0 0 0       = λΦ(x) + t AΦ(x); Φ(x) =       ϕ 1 ϕ 2 φ 1 φ 2 ϕ 0       (41)
with related matrix defined by

A =       -k 0 k 0 0 0 -k 0 k 0 k 0 -k -K 1 0 K 1 0 k 0 -k -g 0 0 0 K 1 g -K 1      
.

Multiplying (40) on the left by t Φ, we deduce

ϕ 1 ∂ x U 1 -ϕ 2 ∂ x U 2 = λ t Φ U + t ΦA U.
Taking the transpose of (41) and multiplying on the right by U, we also have

-∂ x ϕ 1 U 1 + ∂ x ϕ 2 U 2 = λ t Φ U + t ΦA U.
As a consequence, we deduce the relation

(U 1 ϕ 1 ) -(U 2 ϕ 2 ) = 0, ∀x ∈ [0, L]. (42)
Since U 1 (L) = U 2 (L) in [START_REF] Perthame | Transport equations in biology[END_REF] and by initial conditions U 1 (0) = 0, ϕ 2 (0) = 0, then also ϕ 1 (L) = ϕ 2 (L), as set in [START_REF] Shultz | Basic Principles Of Membrane Transport[END_REF]. It means that (U 1 ϕ 1 ) = (U 2 ϕ 2 ) ∀x ∈ [0, L]. Thanks to this relation, we can consider in our previous computation:

ϕ 2 (x) = U 1 (x) U 2 (x) ϕ 1 (x), ∀x ∈ [0, L].

  (a) Concentration profiles with permeability Pi = 2 • 10 -7 [m/s]. (b) Concentrations in 2D with Pi = 2 • 10 -7 . Length of lumen on vertical axis. (c) Concentration profiles with permeability Pi = 2 • 10 -5 [m/s]. (d) Concentrations in 2D with Pi = 2 • 10 -5 . Length of lumen on vertical axis.

Figure 2 :

 2 Figure 2: Concentration profiles for V m,2 = 2πr 2,e 10 -5 and different permeability values.

  (a) Axial concentrations in the lumen for different values of permeability. (b) Fractional increase in concentration as a function of permeability.

Figure 3 :

 3 Figure 3: Concentration profiles for V m,2 = 2πr 2,e • 10 -5 [mol.m -1 .s -1 ].

  (a) Axial concentrations in the lumen for different values of permeability. (b) Percentage of concentration gradient as a function of permeability.

Figure 4 :

 4 Figure 4: Concentration profiles for V m,2 = 2πr 2,e • 10 -4 [mol.m -1 .s -1 ]

Figure 5 :

 5 Figure 5: Percentage of concentration gradient 2D with range V m,2 ∈ 2πr 2,e • (10 -5 , 10 -4 ) [mol.m -1 .s -1 ] (x-axis) and P ∈ 2 • (10 -8 , 10 -5 ) (y-axis)

A Large permeability asymptotic

In this section we consider the case where the permeability between the lumen and the epithelium is large, i.e. when P i → ∞, with i = 1, 2 in the definition of constants k 1 and k 2 . For this purpose, we set

ε and we let ε go to 0. Physically, this means fusing the epithelial layer with the lumen.

Rewriting [START_REF] Feher | Quantitative Human Physiology[END_REF] in this perspective gives

We expect the concentrations u ε 1 and q ε 1 to converge to the same quantity. The same happens for u ε 2 → u 2 and q ε 2 → u 2 . We denote u 1 , respectively u 2 , the limit of u ε 1 and q ε 1 , respectively u ε 2 and q ε 2 . Adding (35a) to (35c) and (35b) to (35d), we obtain

Passing formally to the limit ε → 0, we arrive at

coupled to the equation for the concentration in the interstitium obtained by passing into the limit in equation (35e)

The equations ( 36), (37), (38) describe the same concentration dynamics in a system without epithelium, previously studied in [START_REF] Tournus | Modèles d'échanges ioniques dans le rein: théorie, analyse asymptotique et applications numériques[END_REF] and [START_REF] Tournus | Analysis of a simplified model of the urine concentration mechanism[END_REF]. The formal computation above shows that this 3 × 3 system may be considered as a good approximation of the larger system (6) for large permeabilities. Such a convergence result may be proved rigorously and it is investigated in [15]. It relies on specific a priori estimates and the introduction of an initial layer.