Safaâ Dafrallah

Aouatif Amine

Stéphane Mousset

Abdelaziz Bensrhair

Will Capsule Networks overcome Convolutional Neural Networks on Pedestrian Walking Direction ?*

Thousands of people are dying every year due to road accidents; in fact 23% of world fatal accidents are pedestrians related, where 40% of them occur in Africa as reported by the World Health Organisation (WHO). Predicting the walking direction of a pedestrian could help to avoid an eventual accident. Existing studies can not handle pose and orientation transformations of the input object contrary to our proposed method. This paper describes a novel approach to determine the pedestrian orientation using Capsule Networks (CapsNet) based scheme. CapsNet are a new deep learning architecture that overcome some limitations of the existing studies, they are group of neurons invariant to rotation and affine transformations, which represent a specific interest to this work. Capsule Networks predicts the walking directions of pedestrians to prevent such mortal accidents, using four main walking directions (front, back, left and right).For this purpose, a new pedestrians dataset gathered from the most popular cities in Morocco is collected to be studied and used as a proof of the proposed approach. To enhance this proposed approach, we evaluated it using Daimler dataset and compared it to Convolutional Neural Networks (CNN) architectures.

Experimental results reveal that the performance of the proposed approach reaches an accuracy of 97.60% on daimler dataset and 73.64% on our Moroccan collected dataset.

I. INTRODUCTION

Road traffic injuries are the first cause of death for those who are under 30 years. Approximately 1.35 million person dies annually around the world from road traffic accidents, where 23% of victims are pedestrians [START_REF] Who | Global status report on road safety 2018[END_REF]. According to the World Health Organisation (WHO), 93% of the world fatalities on the road occur in low and middle income countries [START_REF]Road Traffic Injuries[END_REF]. These countries have approximately 60% of the world's vehicles and a high population size.

In order to reduce fatal accidents involving pedestrians, multiple researches are implemented in Advanced Driver Assistance Systems (ADAS) specifically for Pedestrian Crash Avoidance Mitigation (PCAM) [START_REF] Gandhi | Pedestrian collision avoidance systems: A survey of computer vision based recent studies[END_REF]- [START_REF] Chen | Pedestrian-vehicular collision avoidance based on vision system[END_REF]. To that end, researches in PCAM systems focus more on detecting pedestrians than on predicting their walking direction, hence only few ones intend to include pedestrian orientation to those systems. Additionally, existing PCAM systems are designed for well structured areas containing road signs and floor markings, however in low and middle income countries roads are usually poorly structured.

In this paper, Morocco is chosen as a case study, where 28% of road fatal injuries are pedestrian related with 996 deaths in 2016, as reported by the Moroccan Ministry of Environment, Transport, Logistics and Water (METLE) [START_REF] Metle | Recueil des Statistiques des Accidents Corporels de la Circulation de l'année 2016[END_REF].

In this context, we collected a new dataset of moroccan pedestrians from two cities (Rabat and Kenitra). After analyzing this data, we find that traffic laws are less respected by both drivers and pedestrians, where the common way of pedestrians crossing was arbitrarily, specially in poorly structured areas. As a result of this behavior, a high rate of pedestrian accidents is noticed in Morocco.

As a solution, we propose to include pedestrian walking direction to PCAM systems using the moroccan collected dataset. Note that this work is supported by the METLE, in collaboration with the National Center for the Scientific and Technical Research (CNRST).

The remainder of this paper is organized as follows: Section II discussed the Related Work. Section III describes the proposed approach and the collected dataset. Experimental results and discussion are in section IV, and finally the conclusion and future works in section V.

II. RELATED WORK

To predict the pedestrian walking direction the system has first to detect the pedestrian. Starting from hand crafted features to deep learning methods the pedestrian detection field has known remarkable improvements in the last decade.

Dalal et al. [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] proposed a human detection approach based on Histogram of Oriented Gradients (HOG) descriptor, that gives a performance up to 89%. Whereas, Viola et al. [START_REF] Viola | Detecting pedestrians using patterns of motion and appearance[END_REF] used Haar-like Wavelets to detect pedestrians with a rate of 80% and a false positive rate of 1 for every 2 frames. Dollár et al. [START_REF] Dollár | Fast feature pyramids for object detection[END_REF] proposed the Aggregated Channel Features (ACF), that aims to compute 10 channels from an input image, those channels represent the normalized gradient magnitude, Histogram of Oriented Gradients and LUV color. The average miss rate of ACF is about 41%.

Nowadays, automated feature extraction also known by deep learning methods, have shown competitive results for pedestrian detection. As an instance, Bunel et al. [START_REF] Bunel | Detection of pedestrians at far distance[END_REF] used the Convolutional Neural Network (CNN) to detect small scale pedestrians that are at far distance from the camera (30 pixels or less). This method reachs a miss rate lower than 10%. Faster Region-based Convolutional Neural Network (Faster R-CNN) was used by [START_REF] Zhang | Is faster r-cnn doing well for pedestrian detection?[END_REF] for pedestrian detection and achieves an accuracy of 92.7%. While, Redmon et al. [START_REF] Redmon | You only look once: Unified, real-time object detection[END_REF] proposed a new approach for real-time object detection called You Only Look Once(YOLO), which gives an accuracy of 57.9% on COCO dataset but still 1000x faster than R-CNN and 100x faster than Fast R-CNN. Tian et al. proposed in [START_REF] Tian | Pedestrian detection aided by deep learning semantic tasks[END_REF] a multi-task deep model for pedestrian detection that achieves a log-average miss rate of 34.99%.

Once the detection is done, bounding boxes are extracted and used as an input to the pedestrian orientation prediction network. In this point, Gandhi et al. [START_REF] Gandhi | Image based estimation of pedestrian orientation for improving path prediction[END_REF] proposed a Hidden Markov Model (HMM) to model transitions between pedestrian orientations over time and applied the Support Vector Machine (SVM) to estimate the discrete probability distribution of the orientation using 8 orientation bins. This method achives an accuracy of 49.7% for images that fall in the same bins as the ground truth and 81.3% for images that fall in both same or adjacent bins. Shimizu et al. [START_REF] Shimizu | Direction estimation of pedestrian from multiple still images[END_REF] proposed as well an SVM based approach to classify the pedestrians walking directions into 16 directions with a difference of 22.5 degrees. This approach classifies more than 90% of test images into the correct direction. The use of Deep Convolutional Neural Network (DCNN) by Hara et al. [START_REF] Hara | Designing deep convolutional neural networks for continuous object orientation estimation[END_REF] gives an accuracy of 70.6% for a predicted orientation within 22.5 degrees from the ground truth and of 86.1% within 45 degrees. In [START_REF] Dominguez-Sanchez | Recognizing pedestrian direction using convolutional neural networks[END_REF], Sanchez et al. used various state of the art CNN networks such as AlexNet, GoogleNet and ResNet for pedestrian movement direction recognition. This approach first compute the optical flow of the input image to estimate the direction of the pedestrian as one of the three predefined directions (left, right and front), and use afterwards the output image as an input to the CNN architecture. As a result, ResNet achieves the best accuracy of 94% on the validation set.

In spite of the good performing results in pedestrian detection and orientation prediction, Convolutional Neural Network has as well some limitations. In fact, the use of the pooling layer makes CNN losing valuable information about the relative position and the orientation of the object, which makes CNN invariant for object translation and rotation as reported in [START_REF] Lecun | Deep learning[END_REF] by Lecun et al. As a solution, Capsule Networks was proposed by the team of G. Hinton in [START_REF] Sabour | Dynamic routing between capsules[END_REF] achieving a test error of 0.25% on MNIST dataset, they were applied in several research fields such as Traffic Signs detection [START_REF] Kumar | Novel deep learning model for traffic sign detection using capsule networks[END_REF] with an accuracy of 97.6%. To the best of our knowledge, the CapsNet was never used on pedestrian orientation classification.

In this paper, we proposed a novel pedestrian walking direction prediction approach, based on Capsule Networks that takes as an input pedestrian bounding boxes detected using YOLO algorithm since it produces fast and good results on pedestrians detection.

III. PROPOSED APPROACH

The aim of this work is to decrease pedestrian fatal accidents, by making a new ADAS system taking into consideration pedestrians orientation. This system detects pedestrians when moving around the danger zone of the vehicle in poorly structured areas, and predict thereafter the pedestrian direction. In this paper, the effective pedestrian movement has not a proper interest due to moving vehicles, which may give a stable status to the pedestrian and makes his orientation intention more important. The danger zone is defined as an area where the vehicle cannot stop in time. In other words, it is the area that is less than or equal to the stopping distance of the vehicle, where:

Stopping Distance = Reaction Distance + Braking Distance (1)
However, the focus of this paper is on the outer zone, which will be designated as alert zone (see Fig. 1).

According to Massachusetts driver's manual, the twosecond rule is a minimum safe distance for good road conditions and moderate traffic. Whereas, for more safety three to four seconds must be added. Thus, the alert zone of four seconds is used for the system which gives a distance of 64m for a vehicle speed of 60km/h.

An accident between a vehicle and a pedestrian is more likely to occur when the pedestrian walks towards the vehicle within the alert zone. The general scheme of a vehicle and a pedestrian in a danger and alert zone is defined based on the vehicle speed as well as the distance between the vehicle and the pedestrian, which is outlined by the lateral distance called L and the longitudinal distance named R as shown in Fig. 1.

In this study, the system is considered as three major factors: Firstly, the number of people in the nearby zone is considered to be one of the main attributes of the road. Therefore, the areas are going to be divided in high an low density areas. Another major factor is the direction of the vehicle itself, it must deduce the surrounding danger regarding the actual direction of the vehicle whether it is going straight or taking a turn. Finally, the orientation of the pedestrian is regarded as the third major factor. In fact, this work shrinks the pedestrians' moving directions to four major ones, that are moving forward and backward inside the pavement, and crossing the street left to right and right to left. Actually, this first part of the approach mainly emphasises the movement of pedestrians in all the described directions within an environment of low density area and up-front moving vehicle.

This system regards pedestrians walking either within the sidewalks or crossing the road. The first case is considered to be safe; therefore, no intervention is needed. On the other hand, pedestrians crossing the road are more likely to get involved in mortal accidents. Hence, the system should alert the driver to pay more attention to the upcoming danger. Thus, not only the detection of pedestrians is a major requirement, but also the walking direction, which is the main concentration of this approach, is fundamental for this system to work properly. In order to select only pedestrians in risky situations, a pedestrian orientation prediction algorithm based on Capsule Networks is proposed.

In this current section, we will firstly introduce Capsule Networks, then we will describe our moroccan collected dataset used as an input to the network, and finally we will present our proposed architecture for pedestrians' walking direction.

A. Capsule Networks

A capsule is a group of neurons whose activity vectors represent the pose parameters of an entity, and the vectors length represents the existence probability of that entity. Unlike the convolutional network, capsules conserve detailed information about the location and the pose of the entity. A slight rotation of the image involved a slight change in the activation vector. The capsule network as represented in [START_REF] Sabour | Dynamic routing between capsules[END_REF], contains two convolutional layers and one fully connected layer. The first convolutional layer extracts feature maps from the input image. The result is then resized to an array of vectors with components, which constitute the input to the second convolutional layer called primary capsule layer. The vector length must be between 0 and 1 since it represents the existence probability of an entity. Hence, a squash function is applied to shrunk the long vectors to nearly one and the short vectors to almost zero.

The main role of the third layer, which is the fully connected layer, is classification. This latter contains a capsule per class, where each capsule has a belonging probability for all classes and the class having the highest probability, is assigned to the capsule. Capsules from the primary capsule layer (second layer) predict the output vectors of the fully connected layer (third layer). The prediction ûj/i is produced by multiplying the output vector u i of a capsule in the second layer with a transformation matrix W ij :

ûj/i = W ij u i (2)
The transformation matrix is learned by the network gradually using backpropagation in the learning process of the primary capsule layer. Thereafter the agreement a ij between the prediction value made by the capsule i of the second layer ûj/i and the current output vector of the capsule j in the third layer v j , is calculated using a dot product as represented in( 3)

a ij = ûj/i .v j (3)
For each predicted vector, a routing weight is used called b ij and initialized by zero for all capsules in the both layers. Then a softmax function c ij is applied to that routing weight for each capsule in the primary capsule layer. The weighted sum s j of all prediction vectors, is after that calculated for each capsule in the fully connected layer:

s j = i c ij ûj/i (4)
Then a squash function is applied to that weighted sum, giving as a result the real outputs v j of capsules of the third layer:

v j = ||s j || 2 1 + ||s j || 2 s j ||s j || (5)
Subsequently, the routing weight b ij is updated by adding to it the agreement between the predicted and the real vector.

b ij = b ij + ûj/i .v j (6)
This process represents one iteration of the routing algorithm. In case of a right prediction, the routing weight increases which increased the length of the output vector in the next iteration; thus, the existence probability of the entity represented by the vector. According to [START_REF] Sabour | Dynamic routing between capsules[END_REF], the recommended number of routing iterations is three.

The length of the output vector is subsequently used to measure the probability that the entity exists by calculating the margin loss L k . A separate margin loss is used for each class k as shown in [START_REF] Metle | Recueil des Statistiques des Accidents Corporels de la Circulation de l'année 2016[END_REF].

L k = T k max(0, m + -||v k ||) 2 + λ(1 -T k )max(0, ||V k || -m -) 2 (7)
Taking into account that:

• T k = 1 if the entity is present.
• m + and m -are the hyperparameters which equal respectively 0.9 and 0.1. • λ = 0.5.

B. Used datasets 1) Daimler dataset:

We first train the network on Daimler Monocular Pedestrian detection images of 48x96px that we downsampled to 48x48 for computational resource purpose. 12000 samples are used for the training step and 1000 for testing. Daimler dataset doesn't provide ground truth information about the four pedestrian walking directions (front, back, left and right). To fit our purpose, we manually annotated the samples to the four pedestrian orientation classes described above. Despite the high image quality and the neat dataset that daimler provides, the samples are taken from well structured areas of a european city which represents the perfect case for a classification system. To evaluate the proposed architecture on a less structured area, we collected our own dataset gathered from Moroccan cities.

2) Moroccan collected dataset: The dataset contains images obtained from natural scenes of pedestrians walking and crossing the road, gathered from two Moroccan cities, Rabat and Kenitra. The acquisition was done in well structured avenues containing floor marking and pavements, in addition to poor structured ones where the place is crowded and pedestrians walk arbitrarily as seen in Fig. 2. To have a variant database the collection was during different lighting and weather conditions. In poor structured areas and due to the absence of road signs, pedestrians act randomly and cross between moving vehicles, specially in a crowded place. Therefore, in well structured areas known by a fast traffic flow, road signs are less respected by road users. This results in more serious accidents endangering pedestrians lifes.

Our collected data contain approximately 2580 cropped pedestrian images, that are extracted from a three hours video recording of one minute each, captured by a monochrome industrial camera with CMOS sensor and a resolution of 2.3 MP for a maximum of 60 fps. A picture of the camera on board the testing vehicle is shown in Fig. 2. To augment our dataset, images are mirrored using horizontal flip, which doubled the size of the dataset to 5160 pedestrian samples, that are resized afterwards to 48x48 pixels. The dataset contains four classes, each one represents one of the four orientations (front, back, left and right) of the pedestrian walking direction as represented in Fig. 3. For training, 4160 images among the 5160 images are used, where each class contains 1040 samples. While, the testing was done using 250 samples for each class, so a total of 1000 samples for test.

C. Capsnet Architecture for Pedestrian Walking Direction

We build a pedestrian orientation classification system based on capsule networks. The system contains the encoder and decoder part as illustrated in Fig. 4. While the encoder part classifies the input image into one of the four classes, the decoder part reconstructs the input image basing on the result of the classification.

The encoder part contains four layers:

• The two first ones are convolutional layers with 64x5x5 filters of stride 1 for the first one, it takes as an input a 48x48 grayscale image which gives an output of 44x44x64 tensor. While the second one contains 128x5x5 filters with stride of 1 and outputs a tensor of 40x40x128. • The third layer represents the primary capsule layer, it contains 16 channels of 8 dimensions, thus each capsule receives as an input features extracted from the first layer, so a total input of 40x40x128x16. In what concerns the decoder part, this latter contains 3 fully connected layers of 512, 1024 and 2304 filters respectively, it uses the true label from the PedCaps layer to reconstruct the input image.

IV. RESULTS AND DISCUSSION

The proposed approach is tested using both Daimler and our collected dataset giving an accuracy of 97.60% and 73.64% respectively. The gap of the accuracy rate between the two datasets can be merely interpreted by the difference of the pedestrian crossing way for the both datasets. In well structured areas represented by Daimler dataset, pedestrians cross horizontally which makes the classification much easier than in our collected dataset where pedestrians cross obliquely.

We tested the network using different capsule architectures as represented in Table I. To train the network we used 50 epochs of 64 batch sizes. According to the results in Table I, reducing the number of feature maps from 256 to 128 and the number of capsule channels from 32 to 16 using 2 routing iterations, leads to a higher accuracy from 90.62% to 96.87% and a loss of 0.016. However, we choose the last architecture having the best accuracy of 97.60% for Daimler dataset by using two convolutional layers with 64 filters in the first layer and 128 in the second one with 16 primary capsules. Figure 5 illustrates the loss obtained while training the last architecture. The total loss of the algorithm is calculated using the margin loss and the reconstruction loss, with a regularization scale [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF].

T otal loss = M argin loss + λ(Reconstruction loss) [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] With λ = 0.0005 represent the regularization scale. The reconstruction loss is calculated using mean squared error between the input and the reconstructed image. Samples of the reconstructed images are shown in Fig. 6.

We compare the proposed approach to some well known CNN architectures. In Table IV, we present the confusion matrix of the capsnet architecture in Daimler dataset and we compared it with the one published in [START_REF] Dominguez-Sanchez | Recognizing pedestrian direction using convolutional neural networks[END_REF] (see Table V) using 3 orientation classes only. According to the results we can deduct that the use of 4 orientations gives more precision as well as a good classification accuracy comparing to 3 orientations.

Comparing between Daimler and our collected dataset using Capsule Network architecture, we can easily observe

Fig. 1 :

 1 Fig. 1: General scheme of a vehicle and a pedestrian in a danger and alert zone

Fig. 2 :Fig. 3 :

 23 Fig. 2: sample of acquisitions from Kenitra and Rabat captured using a camera on board a moving vehicle

Fig. 4 :Fig. 5 :

 45 Fig. 4: Pedestrian Capsnet Architecture containing the encoder part that classify the input image into one of the four walking directions, and the decoder part that reconstructs the input image using the result of the encoder part

Fig. 6 :

 6 Fig. 6: Samples of reconstructed images

•

  The final layer which we named PedCaps (inspired from capsules for pedestrians) consists of 4 capsules of 16 dimension, each capsule refers to one class among the four orientation classes. This layer aims to classify the input image and assign it to one of the classes aforementioned.

  Table II shows the test accuracy of Alexnet and Resnet on Daimler dataset that achieved an accuracy of 95.52% and 96.45% respectively, while the

  Table III represents confusion matrices of the architectures presented above on Daimler dataset, we can observe that the best performer network is Capsule networks.

	TABLE II
	ALEXNET, RESNET AND CAPSNET ACCURACY
	COMPARAISON ON DAIMLER DATASET
	Architecture	Accuracy
	Alexnet	95.52%
	Resnet	96.45%
	Capsnet	97.60%
	proposed Capsule Network architecture performed the best
	result of 97.60%.

ACKNOWLEDGMENT

This research work is supported by METLE and CNRST under title: "SAFEROAD Meta-plateforme pour la Sécurité Routière (MSR)" project with contract No. 24/2017, and by LITIS laboratory in INSA Rouen.

As this approach is intended to be involved in a pedestrian collision prediction system, it has to be suitable for real life. Thereby, the use of more than 4 orientation bins is highly recommended, and will be the intent of future works.

V. CONCLUSIONS

This paper tackles one of the most crucial issues, which is pedestrians' walking directions, by presenting a novel approach based on Capsule Networks. This approach aims to classify pedestrians' orientations into one of the main four directions (front, back, left and right), the proposed CapsNet architecture performed better results for pedestrian orientation classification with an accuracy of 97.60% on Daimler dataset, compared to the tested CNN architectures. Additionally, we proposed in this paper a new dataset captured from Moroccan cities using a fixed camera on-board a moving vehicle. For future work, we aim to predict the pedistrian direction in video sequences to be applied in the Moroccan area. The integration of this method on PCAM systems is also planned to estimate collision probability and generate appropriate warnings to the Moroccan driver.


F Smaï

A thermodynamic formulation for multiphase compositional flows in porous media

Keywords: porous medium, multiphase flow, thermodynamics, capillarity, interface

The main purpose of this work is to develop a practical and consistent model of multiphase flow in porous media. The thermodynamics of irreversible processes is applied to the representative elementary volume of a porous medium filled with multiple compositional fluids. The entropy production is derived and analyzed in order to establish phenomenological equations for the macroscopic fluxes and an equation of state for the matrix-fluids-interfaces system at equilibrium. The entropy dissipation due to the transformation of the interfaces is identified and its analysis gives a new model for capillarity with two phases or more. The physical model is formulated in a well-posed system of elliptic-parabolic differential equations based on a set of persistent primary variables. The formulation encompasses a wide range of classical models and extends them by handling an arbitrary number of phases and components.

Introduction

The idea that thermodynamics can help in building macroscopic descriptions of flows in porous media is explored since several decades [START_REF] Coussy | Mechanics and physics of porous solids[END_REF][START_REF] Everett | Thermodynamics of multiphase fluids in porous media[END_REF][START_REF] Hassanizadeh | Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries[END_REF][START_REF] Hassanizadeh | Thermodynamic basis of capillary pressure in porous media[END_REF][START_REF] Marle | On macroscopic equations governing multiphase flow with diffusion and chemical reactions in porous media[END_REF][START_REF] Morrow | Physics and thermodynamics of capillary action in porous media[END_REF]. It led to significant insights, like interpreting the empirical Darcy's law [START_REF] Darcy | Les Fontaines publiques de la ville de Dijon[END_REF] as the first order approximation of a dissipative process [START_REF] Marle | On macroscopic equations governing multiphase flow with diffusion and chemical reactions in porous media[END_REF]. The purpose of thermodynamic approaches is to give theoretical foundations to the most used macroscopic equations and to develop accurate models for less known processes. One of their strengths is consistency, in opposition to a modelling strategy based on the analogy with phenomena in continuous medium and the superposition of the corresponding equations.

Although the thermodynamic approaches are multiple, they have in common to start with a thermodynamic description at the pore level, then to derive a macroscopic description from them. The approaches can differ at both steps. Concerning the transition from pore level to macroscopic level, the challenge is to properly relate the quantities and relations used at the two levels. The notion of a length scale separating between two levels of description appears independently in porous media and continuum thermodynamics. For porous media, it corresponds to the representative elementary volume [START_REF] Bear | Dynamics of fluids in porous media[END_REF] and it separates between coarse and smooth descriptions of the porous medium. For continuum thermodynamics, it corresponds to the elemental volume on which the local equilibrium is valid [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF][START_REF] Kondepudi | Modern thermodynamics: from heat engines to dissipative structures[END_REF] and it separates between uniform and non-uniform description of the continuum. These two volumes both have their own length scale that depends on the structure of porous medium or the processes occurring inside. The hierarchy between the two length scales is not known a priori, they can be the same or one can be larger than the other. Nevertheless applying thermodynamic principles to porous media involves both length scales and the way they compare has a direct impact on the final macroscopic description.

One goal of modelling flows in porous media is the ability to perform numerical simulations, whose importance is constantly increasing in earth sciences and engineering. The expected advances concern the speed and accuracy of the simulations [START_REF] Chen | Computational methods for multiphase flows in porous media[END_REF] but also their complexity by coupling multiple phenomena [START_REF] Kolditz | Thermo-hydro-mechanicalchemical processes in porous media[END_REF]. While improving the numerical methods is an active topic, the models used for simulating flows in porous media are generally not based on thermodynamics: Darcy's law is initially an empirical law and its extensions are mainly based on hydrodynamic arguments; the relative permeability that adapts Darcy's law for multiphase flow lacks theoretical foundations and satisfactory agreement with experiment; heat transfer, diffusion or chemical reactions are generally superimposed to the flow, using equations obtained by analogy with those governing the same phenomena in a continuous medium. Thus these models are exposed to some thermodynamic inconsistencies that can be the root of inherent biases and numerical difficulties. Concerning the coupling with other phenomena, the thermodynamic modelling has proven its relevance for the interactions between poro-mechanics and flows [START_REF] Coussy | Mechanics and physics of porous solids[END_REF]. Reactive transport is an other active topic that could exploit a proper modelling of compositional flows [START_REF] Sin | Multiphase multicomponent reactive transport and flow modeling[END_REF]. Finally, numerous formulations [START_REF] Abadpour | Method of negative saturations for modeling two-phase compositional flow with oversaturated zones[END_REF][START_REF] Bourgeat | Two-phase, partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository[END_REF][START_REF] Bourgeat | On persistent primary variables for numerical modeling of gas migration in a nuclear waste repository[END_REF][START_REF] Class | Numerical simulation of non-isothermal multiphase multicomponent processes in porous media.: 1. an efficient solution technique[END_REF][START_REF] Coats | An equation of state compositional model[END_REF][START_REF] Forsyth | A two-phase, two-component model for natural convection in a porous medium[END_REF][START_REF] Jaffré | Henry'law and gas phase disappearance[END_REF][START_REF] Neumann | Modeling and simulation of twophase two-component flow with disappearing nonwetting phase[END_REF] have been developed to address the particularities of multiphase flow, like appearance and disappearance of fluid phases or switching between parabolic, elliptic and hyperbolic behaviors. They collectively show how to deal with multiphase flows but they are ad hoc constructions, difficult to generalize. A mathematical formulation of flows in porous media derived from thermodynamic principles would be more generic, ensure consistency, hopefully also when properly coupled with other phenomena, and could moreover benefit from the mathematical tools developed over times for thermodynamic systems.

This work proposes a derivation of macroscopic equations for multiphase flow in porous media. A strong assumption is that local equilibrium applies at the representative elementary volume. This means there is a length scale at which the thermodynamic quantities are uniform and the porous medium is seen as a smooth continuum. Although it is mostly suited to slow flows that are close to equilibrium, this assumption proves to be sufficient to recover most of the classical models. At first, the porous system is divided into distinct parts and the classical relations of equilibrium thermodynamics are reminded for each one. Next, the balance equation for entropy is derived from the matter and energy conservation principles and the local equilibrium of each part of the porous system. The entropy production is then analyzed to identify the different processes and, with additional assumptions, phenomenological laws are formulated for each processes. The resulting model describes an arbitrary number of phases, made of an arbitrary number of components, flowing through the porous media, subject to dissolution and to changes of phase, and constrained by capillarity. The analysis of the entropy production additionally provides new definitions for the equilibrium of the system of interfaces and for the pressure-saturation relation that hold for systems with two phase or more. These definitions avoid involving ambiguous quantities like average interface curvature or average contact angle and they introduce new constraints for the pressure-saturation relation. In order to demonstrate the practicality of the full model, a mathematical formulation of the model is developed. It takes the form of a system of elliptic-parabolic differential equations whose main properties are: a persistent set of primary variables, a symmetric positive structure and a trivially associated entropy equation.

Description of the porous medium system

This work adopts a description of the porous medium and its content as a continuum. At each point of this macroscale continuum, the different parts of the porous system are independently at local equilibrium and are governed by their own equations of state. This section details the assumptions and the notations of that description.

Partitioning the system

Any volume, V, of the porous medium system is assumed to be divided between the volume of the rock matrix (the solid phase), V r , and the volume of each fluid phase α ∈ P, V α . Here, P denotes the set of indices for all the possible fluid phases. The fluid phase α ∈ P is said missing in the system volume V when V α = ∅. In the following, the fluids are assumed to always completely fill the pores: V = k∈P∪{r} V k . The measure of these volumes are denoted by V = |V| and V k = |V k | for k ∈ P ∪ {r}. The interface area between two fluid or solid phases indexed by k and l is denoted by A {k,l} = |∂V k ∩ ∂V l |. In the following, I(E) = {{k, l} | k ∈ E, l ∈ E, k l} will denote the set of all pairs of the given index set E; the index set of all interfaces is I(P ∪ {r}). For a fixed system volume V, the volume V α of the fluid α ∈ P can change as the fluids flow through the porous medium. However, the solid phase will be assumed rigid in this work, meaning that V r can not change when V is fixed. According to these definitions, the interfaces have no volume and

V = V r + α∈P V α .
Each phase is made of one or more components and the same component can be found in different phases. C denotes the set of indices for all the possible components and N i denotes the total number of mole of the component i ∈ C in the system volume V. In this work, all the components are assumed to be exclusively in the fluid phases. This means that components in the interfaces are neglected and that there is no exchange of matter with the rock matrix. By denoting N i α the number of mole of component i ∈ C in the volume V α of fluid phase α ∈ P, one get the relation

N i = α∈P N i α ∀ i ∈ C .
Let U denote the internal energy contained in the system volume V. The internal energy is defined as the total energy excluding the kinetic energy due to the motion of the whole system and the potential energy due to external force fields [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF]. As such, the internal energy of the system is the sum of the internal energy of its parts,

U = k∈P∪{r} U k + {k,l}∈I(P∪{r}) U {k,l} ,
where U k denotes the internal energy in the volume V k of the phase k ∈ P∪{r} and U {k,l} denotes the internal energy in the surface ∂V k ∩ ∂V l of the interface {k, l} ∈ I(P ∪ {r}).

Equilibrium thermodynamics

For a given isolated thermodynamic system, the state of thermodynamic equilibrium is the time-invariant state in which we see no further physical or chemical change in the system [START_REF] Kondepudi | Modern thermodynamics: from heat engines to dissipative structures[END_REF]. Equilibrium thermodynamics describe the changes of state of the system between successive equilibria. In particular, the state of the system is specified in terms of state variables such as volume, mole numbers of components, internal energy or entropy. The first law of thermodynamics, that expresses the principle of energy conservation in term of total differential of functions of state, is the basis of the theory. This section presents the assumptions and the relations that describe the thermodynamic equilibrium of the different parts of the porous system.

Fluid phases

The first law of thermodynamics applied to a fluid phase gives

dU α = T α dS α -p α dV α + i∈C µ i α dN i α ∀ α ∈ P , (1)
where S α denotes the entropy of the fluid phase α in the volume V α . The relation [START_REF] Abadpour | Method of negative saturations for modeling two-phase compositional flow with oversaturated zones[END_REF] indicates that the internal energy of the fluid α is a function of the form U α (S α , V α , (N i α ) i ). It also provides the physical meaning for the derivatives of this function: T α is the temperature of the fluid α; p α is the pressure of the fluid α; and µ i α is the chemical potential of the component i in the fluid α. Similarly to the volume V α , the mole numbers N i α or the internal energy U α , the entropy S α is an extensive variable whereas the temperature T α , the pressure p α and the chemical potentials µ i α are intensive variables. Note that since the fluid in V α is at equilibrium, its intensive variables have uniform values across all the volume V α . By applying the Euler's homogeneous function theorem on relation (1) one get

U α = T α S α -p α V α + i∈C µ i α N i α ∀ α ∈ P . (2)
By taking the differential of equation ( 2) and subtracting equation (1), one derives the Gibbs-Duhem relation

0 = S α dT α -V α dp α + i∈C N i α dµ i α ∀ α ∈ P . (3)
Let Ω α denotes the grand (canonical) potential of the fluid α in the volume V α defined by

Ω α = U α -T α S α - i∈C µ i α N i α ∀ α ∈ P , (4)
from equation (1), its differential is

dΩ α = -S α dT α -p α dV α - i∈C N i α dµ i α ∀ α ∈ P .
The extensive variables can be associated to volume quantities. Let c i α be the molar concentration of component i ∈ C in the fluid phase α ∈ P, u α the volume internal energy of the fluid phase α ∈ P, s α the volume entropy of the fluid phase α ∈ P and ω α the volume grand potential of the fluid phase α ∈ P. These volume quantities are linked to theirs extensive counterparts by

N i α = V α c i α , U α = V α u α , S α = V α s α , Ω α = V α ω α ∀i ∈ C ∀α ∈ P .
Using the volume quantities, the Gibbs-Duhem relation (3) can be written as

dp α = s α dT α + i∈C c i α dµ i α ∀α ∈ P . (5)
The volume of fluid α can be scaled to the volume of the whole system to define the volume fraction φ α by

φ α = V α V ∀α ∈ P .
The amount of internal energy of the fluid α per unit of system volume is

U α /V = φ α u α .
According to equation (2), we have

φ α u α = φ α T α s α -φ α p α + i∈C φ α µ i α c i α ∀α ∈ P ,
the differential of this relation gives

d(φ α u α ) = T α d(φ α s α ) -p α dφ α + i∈C µ i α d(φ α c i α ) + φ α s α dT α -φ α dp α + i∈C φ α c i α dµ i α ,
which, using equation ( 5), reduces to

d(φ α u α ) = T α d(φ α s α ) -p α dφ α + i∈C µ i α d(φ α c i α ) ∀α ∈ P . (6)
Similarly, the amount of grand potential of the fluid α per unit of system volume is

Ω α /V = φ α ω α . According to equation (4), its differential satisfies d(φ α ω α ) = d(φ α u α ) -d(φ α T α s α ) - i∈C d(φ α c i α µ i α ) ,
which, using equation ( 6), reduces to

d(φ α ω α ) = -φ α s α dT α -p α dφ α - i∈C φ α c i α dµ i α ∀α ∈ P . (7)

Interfaces

The interfaces are assumed to only exchange heat and work but not matter. Moreover, considering one kind of interface {k, l} ∈ P ∪ {r} in the system volume V, we assume that the set of the interfaces between phases k and l is at equilibrium as a whole. Then, the first law of thermodynamics applied to the interface gives dU {k,l} = T {k,l} dS {k,l} + γ {k,l} dA {k,l} ∀ {k, l} ∈ I(P ∪ {r}) ,

where S {k,l} denotes the entropy of the interface {k, l} in the volume V. The relation [START_REF] Chen | Computational methods for multiphase flows in porous media[END_REF]

and the Gibbs-Duhem relation for the interfaces is obtained by subtracting the equation [START_REF] Chen | Computational methods for multiphase flows in porous media[END_REF] to the differential of the equation ( 9) 0 = S {k,l} dT {k,l} + A {k,l} dγ {k,l} ∀ {k, l} ∈ I(P ∪ {r}) .

Let Ω {k,l} denotes the grand potential of the interface {k, l} in the system volume V defined by Ω {k,l} = U {k,l} -T {k,l} S {k,l} ∀ {k, l} ∈ I(P ∪ {r}) ,

from equation [START_REF] Chen | Computational methods for multiphase flows in porous media[END_REF], its differential is dΩ {k,l} = -S {k,l} dT {k,l} + γ {k,l} dA {k,l} ∀ {k, l} ∈ I(P ∪ {r}) .

Surface quantities are introduced by scaling the extensive variables to interface area. For each interface {k, l} ∈ I(P ∪ {r}), let u {k,l} be the surface internal energy, s {k,l} the surface entropy and ω {k,l} the the surface grand potential. The surface quantities satisfy U {k,l} = A {k,l} u {k,l} , S {k,l} = A {k,l} s {k,l} , Ω {k,l} = A {k,l} ω {k,l} ∀ {k, l} ∈ I(P ∪ {r}) .

Using the surface quantities, the Gibbs-Duhem relation [START_REF] Class | Numerical simulation of non-isothermal multiphase multicomponent processes in porous media.: 1. an efficient solution technique[END_REF] can be written as dγ {k,l} = -s {k,l} dT {k,l} ∀ {k, l} ∈ I(P ∪ {r}) .

(

) 12
The area of interface {k, l} can be scaled by the volume of the whole system to define the density of interface area a {k,l} as

a {k,l} = A {k,l} V ∀ {k, l} ∈ I(P ∪ {r}) .
The amount of internal energy of the interface {k, l} per unit of system volume is

U {k,l} /V =

Matrix

The rock matrix is assumed to be a non deformable phase, it only exchange heat but neither work nor matter. Under these assumptions, the first law of thermodynamics applied to the rock matrix gives dU r = T r dS r [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF] where S r denotes the entropy of the rock matrix in the system volume V and T r its temperature. A common way to describe a rigid solid is by using the concept of heat capacity. Let c v be the volume heat capacity at constant volume of the rock matrix defined as

c v (T r ) = 1 V ∂U r ∂T r V . ( 16
)
The grand potential of the rock matrix is defined by

Ω r = U r -T r S r
and, from equation ( 15), its differential is

dΩ r = -S r dT r . (17)
Let u r be the volume internal energy of the rock matrix and s r the volume entropy that satisfy U r = V r u r , S r = V r s r .

Besides, let φ r be the volume fraction of rock matrix such that

φ r = V r V .
Using φ r , u r and s r , equations ( 15) and ( 16) give

φ r u r = T r 0 c v (T )dT and φ r s r = T r 0 c v (T ) T dT . ( 18
)
For a fixed system volume V, one can derive from equations ( 15) and ( 17) that

d(φ r u r ) = T r d(φ r s r ) (19)
and d(φ r ω r ) = -φ r s r dT r .

(20)

Representative elementary volume and the continuum approach

Thermodynamics manages to describe systems with a large number of molecules by abandoning the molecular level description in favor of a larger scale description. As a consequence, there is a low limit on the size of a thermodynamic system. The same idea lies behind the concept of representative elementary volume (REV) for porous media [START_REF] Bear | Dynamics of fluids in porous media[END_REF]. At the pore scale, a porous medium exhibits a complex geometry and fluid volumes are difficult to track. A volume V of the porous medium system could be considered a REV if it is smooth enough at the pore scale (e.g. a sphere) and if it is much larger than a single pore. For instance, regarding the porosity, V must be large enough to minimize the fluctuations of the matrix volume fraction due to the random pore size distribution, however V must also be small enough to lower the sensitivity of the volume fraction to large scale inhomogeneity. The same reasoning holds for other averaged quantities such as the content of mole of a component or the content of energy. This work assumes the existence a common concept of REV which remains relevant for all the considered averaged quantities. By introducing the concept of REV, the continuum approach replaces the actual medium and its pore scale heterogeneity by a fictitious continuum in which the averaged quantities have a value at any point of the domain.

Local equilibrium

As heat and fluids flow through the porous medium, the thermodynamic equilibrium can not be assumed for the porous medium system over the entire domain. Yet, thermodynamic quantities such as temperature or pressure remain relevant locally. The local equilibrium assumption states that equilibrium thermodynamic relations are valid for the thermodynamic variables assigned to an elemental volume [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF][START_REF] Kondepudi | Modern thermodynamics: from heat engines to dissipative structures[END_REF]. As a consequence, all intensive thermodynamic variables and volume quantities (densities) becomes functions of position and time. The local equilibrium is an efficient approximation for a wide range of hydrodynamic and chemical systems. The present work adopts this assumption in the context of porous media by using the concept of REV (section 2.3) as elemental volume.

In order to apply both the continuum approach and the local equilibrium to the equilibrium thermodynamics presented in section 2.2, we only retain the equations regarding intensive or averaged variables; namely equations ( 5), ( 6) and [START_REF] Brooks | Hydraulic properties of porous media[END_REF] for the fluids, equations ( 12), ( 13) and ( 14) for the interfaces, and equations ( 18), [START_REF] Everett | A general approach to hysteresis. part 2: Development of the domain theory[END_REF] and [START_REF] Everett | A general approach to hysteresis[END_REF] for the matrix. According to the local equilibrium, all these equations are considered locally valid.

Equations of state of the parts of the system

An equation of state relates the state variables of a thermodynamic system such as its pressure, its temperature or its volume. Equations of state cannot be derived from the general relations of thermodynamics alone but require additional knowledge either from experimental measurements or theoretical considerations on the microscopic structure. An equation of state fully characterizes the thermodynamic behaviour of a substance however it can take various equivalent forms depending on the subset of state variables selected to express it.

Equation ( 18) is an equation of state for the rock matrix, provided that the function c v (T ) is known. For each fluid phase α, the Gibbs-Duhem relation ( 5) suggests an equation of state of the form p α (T α , (µ i α ) i ) . Similarly, equation ( 12) suggests the form

γ {k,l} (T {k,l} )
for the equation of state of the interface {k, l}. The proposed forms of equations of state rely on the assumptions made in this section on each part of the porous medium system, such as the matrix being non deformable or the interfaces not exchanging matter. It should be noted that knowing the equations of states for each part of the system is not sufficient to characterize the whole system. In particular, the volume fractions, φ α , and the densities of interface area, a {k,l} , remain uncontrolled.

Balance equations

The section 2 established a local description of the porous medium system. The space and time evolution of the system will be derived from conservation principles. Following the local description, the balance equations are expressed locally, in the form of partial derivative equations, for each part of the porous medium system.

Conservation of components

The general balance of component i ∈ C in the the fluid phase α ∈ P is

∂φ α c i α ∂t + div ϕ n,i α = Λ n,i α
where ϕ n,i α is the molar flux of component i in the fluid phase α and Λ n,i α is the molar production rate of component i in fluid phase α. The local production of component i in phase α can come from other species of the phase α or from other parts of the system (other fluids, matrix and interfaces). All these transformations can be described as chemical reactions [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF][START_REF] Kondepudi | Modern thermodynamics: from heat engines to dissipative structures[END_REF], the local component production being expressed in terms of reaction velocities v j and stoichiometric coefficients b j,i α :

Λ n,i α = j∈ reactions b j,i α v j .
In the current work, we are not taking in account the full range of possible chemical reactions but only component transfer from fluid phase to fluid phase. The local production terms then find a simpler expression,

Λ n,i α = β∈P Λ n,i α,β
where Λ n,i α,β is the molar rate of component i received by fluid phase α from fluid phase β. By definition of the exchange term Λ n,i α,β , one have

Λ n,i α,β = -Λ n,i β,α ∀ i ∈ C ∀ α ∈ P ∀ β ∈ P . (21)
Furthermore, according to the assumptions of section 2, components are only present in the fluid phases, thus the conservation of components is fully described by the balances of component i ∈ C in fluid phase α ∈ P:

∂φ α c i α ∂t + div ϕ n,i α = β∈P Λ n,i α,β (22)
Let M i denotes the molar mass of the component i ∈ C, and ρ α denotes the mass density of the phase α ∈ P, they satisfy

ρ α = i∈P M i c i α . (23)
We introduce the notion of mass average velocity [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Bird | Introductory transport phenomena[END_REF] of the phase α, q α , defined by

ρ α q α = i∈C M i ϕ n,i α , (24)
where ρ α q α represents the momentum per unit volume of the phase α. The part of the molar flux not due to the phase velocity is defined by

J n,i α = ϕ n,i α -c i α q α . (25)
J n,i α is called the molar diffusive flux (relative to mass average velocity) of the component i ∈ C in the fluid phase α ∈ P. From equations ( 23), ( 24) and ( 25), the diffusive fluxes in the fluid α satisfy

i∈C M i J n,i α = 0 ∀ α ∈ P .
The total balance equation for component i ∈ C is obtained by summing equation ( 22) over all the fluid phases. Using equations ( 21) and (25), we get

∂n i ∂t + div         α∈P c i α q α + J n,i α         = 0 ∀ i ∈ C (26)
with

n i = α∈P φ α c i α ∀ i ∈ C . ( 27
)
The continuity equation is obtained by multiplying the equations ( 26) by the molar mass M i then summing over the components:

∂ ∂t         α∈P φ α ρ α         + div         α∈P ρ α q α         = 0 .

Balance of internal energy

The principle of energy conservation does not apply to internal energy but to total energy (sum of internal energy, kinetic energy and potential energy). Nonetheless a balance equation for the internal energy can be obtained at the microscopic scale by subtracting the balance equation of mechanical energy to the total energy conservation [START_REF] Bird | Introductory transport phenomena[END_REF][START_REF] De Groot | Non-equilibrium thermodynamics[END_REF]. Averaging the internal energy balance over the REV is non trivial, especially because of the non linear contributions of the velocity, and several authors have proposed macroscopic balance equations for the internal energy [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Cheng | Heat transfer in geothermal systems[END_REF][START_REF] Forsyth | Three-dimensional modelling of steam flush for dnapl site remediation[END_REF][START_REF] Lake | Enhanced oil recovery[END_REF][START_REF] Vafai | Handbook of porous media[END_REF]. Unfortunately, the various propositions are not equivalent, neither are they all compatible. Here, we use a simplified balance of internal energy where the only remaining velocity-related contributions are all linear. This is relevant for flow regimes where velocities at pore scale are small enough to neglect high order contributions.

Considering the creeping flow regime [START_REF] Bird | Introductory transport phenomena[END_REF] at the REV scale, we can neglect inertia effects like kinetic energy and viscous dissipation against internal energy, compression effects, heat diffusion and gravity work. According to this assumption and those of section 2, the balance of internal energy for the fluid phase α ∈ P takes the form:

∂φ α u α ∂t + div(h α q α + J u α ) = p ∈P∪{r}∪I(P∪{r}) Λ u α,p + ρ α q α • g . ( 28
)
Where h α is the volume enthalpy of fluid phase α, J u α is the heat flux through fluid phase α, Λ u α,p is the volume rate of energy received by fluid phase α from the part p of the porous system, and g is the gravitational acceleration. The volume enthalpy of the fluid α is related to its internal energy and its pressure by

h α = u α + p α ,
meaning that the term div(h α q α ) of equation ( 28) gathers two contributions: the internal energy advection, div(u α q α ), and the reversible compression work, div(p α q α ). The last term of equation ( 28), ρ α q α • g, is the reversible work of the force of gravity on the fluid α.

As the rock matrix is assumed to be no deformable and to only exchange heat but not work, its balance of internal energy is simpler to establish and take the form

∂φ r u r ∂t + div(J u r ) = p ∈P∪{r}∪I(P∪{r}) Λ u r,p . (29)
Where J u r is the heat flux through the solid phase, and Λ u r,p is the volume rate of energy received by the solid phase from the part p of the porous system.

Since the transport by macroscopic displacement of the interfaces is not considered here, the balance of internal energy for the interface {k, l} ∈ I(P ∪ {r}) take a similar form ∂a {k,l} u {k,l} ∂t

+ div(J u {k,l} ) = p ∈P∪{r}∪I(P∪{r}) Λ u {k,l},p . (30)
Where J u {k,l} is the heat flux through the interface {k, l}, and Λ u {k,l},p is the volume rate of energy received by the interface {k, l} from the part p of the porous system.

In equations ( 28), ( 29) and ( 30), the terms Λ u p,p encompass all the energy exchanges that occur between the parts p and p of the system P ∪ {r} ∪ I(P ∪ {r}). The origin of these exchanges can be heat, work, matter or a combination of these. By definition of the exchange terms, it holds that

Λ u p,p = -Λ u p ,p ∀ p ∈ P ∪ {r} ∪ I(P ∪ {r}) ∀ p ∈ P ∪ {r} ∪ I(P ∪ {r}) . ( 31
)
Thanks to [START_REF] Marle | On macroscopic equations governing multiphase flow with diffusion and chemical reactions in porous media[END_REF], summing the equations ( 28), ( 29) and [START_REF] Leverett | Capillary behavior in porous solids[END_REF] gives the following balance of internal energy for the mixture

∂u ∂t + div          α∈P h α q α + p∈P∪{r}∪I(P∪{r}) J u p          = α∈P ρ α q α • g (32)
with

u = α∈P φ α u α + {k,l}∈I(P∪{r})
a {k,l} u {k,l} + φ r u r .

Balance of entropy

When considering an arbitrary thermodynamic system with entropy S , the variation of entropy dS can be written as the sum of the entropy supplied to the system by its surrounding, d e S , and the entropy produced inside the system, d i S : dS = d e S + d i S . The second law of thermodynamics states that d i S ≥ 0 for any transformation, the case d i S = 0 corresponding to a reversible transformation. Following the continuum approach, this can be expressed locally and the local entropy production by unit volume must also be positive [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF][START_REF] Kondepudi | Modern thermodynamics: from heat engines to dissipative structures[END_REF]. For the porous medium system, the local form of the second law of thermodynamics may be written as

∂s ∂t + div ϕ s = Ψ ≥ 0 ( 34
)
where s is the total volume entropy of the system, ϕ s is the total entropy flux, and Ψ is the rate of local entropy production. The total entropy production Ψ includes the dissipative processes occurring in each part of the system but also those resulting from the interactions between the parts of the system; its expression is derived below.

Entropy being an extensive quantity, the total volume entropy of the system is the sum over the parts of the system:

s = φ r s r + α∈P φ α s α + {k,l}∈I(P∪{r}) a {k,l} s {k,l} . (35)
By taking the time derivative of the sum and using the relations ( 6), ( 13) and ( 19), we get ∂s ∂t = 1 T r ∂φ r u r ∂t

+ α∈P         1 T α ∂φ α u α ∂t + p α T α ∂φ α ∂t - i∈C µ i α T α ∂φ α c i α ∂t         + {k,l}∈I(P∪{r}) 1 T {k,l} ∂a {k,l} u {k,l} ∂t - γ {k,l} T {k,l}
∂a {k,l} ∂t .

By using the notations

       τ p ≡ -1/T p ∀p ∈ P ∪ {r} ∪ I(P ∪ {r}) ν i α ≡ µ i α /T α ∀α ∈ P ∀i ∈ C (36)
and substituting the energy and mole balance equations ( 28), ( 30), ( 29) and ( 22), it comes ∂s ∂t =p,p ∈P∪{r}∪I(P∪{r})

τ p Λ u p,p - i∈C α,β∈P ν i α Λ n,i α,β + p∈P∪{r}∪I(P∪{r}) τ p div J u p + i∈C α∈P ν i α div J n,i α α∈P         τ α div(h α q α ) + i∈I(P∪{r}) ν i α div(c i α q α ) -τ α ρ α g • q α         + α∈P p α T α ∂φ α ∂t - {k,l}∈I(P∪{r})
γ {k,l} T {k,l} ∂a {k,l} ∂t .

The equation can be simplified and written in the form (34) by using the divergence of the scalar-vector product, the anti-symmetry of the Λ • p,p exchange coefficients and the fact that the relation ( 5) is equivalent to

1 T α dp α = h α dτ α + i∈C c i α dν i α .
Finally, the entropy flux is given by

ϕ s = - p∈P∪{r}∪I(P∪{r}) τ p ϕ u p - α∈P i∈C ν i α ϕ n,i α (37)
where

ϕ u α = h α q α + J u α ϕ u {k,l} = J u {k,l} ϕ u r = J u r ,
and the entropy production is given by

Ψ = - 1 2 p,p ∈P∪{r}∪I(P∪{r}) Λ u p,p (τ p -τ p ) - 1 2 i∈C α,β∈P Λ n,i α,β (ν i α -ν i β ) - p∈P∪{r}∪I(P∪{r}) J u p • ∇ τ p - α∈P i∈C J n,i α • ∇ ν i α - α∈P 1 T α q α • (∇ p α -ρ α g) + α∈P p α T α ∂φ α ∂t - {k,l}∈I(P∪{r})
γ {k,l} T {k,l} ∂a {k,l} ∂t .

(

) 38
The expression (37) of ϕ s is similar to the usual form found in general literature [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF][START_REF] Kondepudi | Modern thermodynamics: from heat engines to dissipative structures[END_REF] while the equation [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF] for the entropy production Ψ is very specific to flows in porous media. The next section is devoted to the analysis of the differents terms of Ψ.

Entropy production: equilibria and linear regime

An expression [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF] of the entropy production, Ψ ≥ 0, has been derived in the previous section for a wide range of flows in porous media. The positivity of Ψ constrains the set of acceptable phenomenological laws that can be chosen to model flows. With additional assumptions, analysing the entropy production provides more precise constrains. This section studies the consequences of two simple modelling ideas: (i) the flow is slow enough to assume the different parts of the system are locally at equilibrium; (ii) the flow is slow enough to assume the system is spatially close to equilibrium and follows the linear regime.

The symmetry principle

The symmetry principle [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF][START_REF] Kondepudi | Modern thermodynamics: from heat engines to dissipative structures[END_REF] is helpful to split the entropy production in independent terms. It is based on the idea that the causes cannot have more symmetries than their effects. Applied to the entropy production, it suggests to group the terms of the sum according to their symmetries and to impose the positivity for each group. For equation [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF], three kinds of term are involved, each with a different kind and amount of symmetries: the terms with a scalar difference, the terms with a gradient and the terms with a time derivative. Rewriting the entropy production with these three groups gives the following decomposition:

Ψ = Ψ difference + Ψ gradient + Ψ time derivative Ψ difference = - 1 2 p,p ∈P∪{r}∪I(P∪{r}) Λ u p,p (τ p -τ p ) - 1 2 i∈C α,β∈P Λ n,i α,β (ν i α -ν i β ) ≥ 0 Ψ gradient = - p∈P∪{r}∪I(P∪{r}) J u p • ∇ τ p - α∈P i∈C J n,i α • ∇ ν i α - α∈P 1 T α q α • (∇ p α -ρ α g) ≥ 0 Ψ time derivative = α∈P p α T α ∂φ α ∂t - {k,l}∈I(P∪{r}) γ {k,l} T {k,l} ∂a {k,l} ∂t ≥ 0 . (39)

The linear regime

The entropy production terms Ψ difference and Ψ gradient of equation ( 39) have the general form

Ψ x = y J x y F x y ≥ 0
where J x y and F x yz are respectively called flows and forces in a general thermodynamics sense. At equilibrium, all the forces and the corresponding flows vanish. For a system close to equilibrium, the flows can be expected to be linear functions of the forces [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF][START_REF] Kondepudi | Modern thermodynamics: from heat engines to dissipative structures[END_REF]. Accordingly, the following relation between the flows and the forces is assumed:

J x y = z L x yz F x z .
Here the coefficients L x yz are called phenomenological coefficients and the matrix they form, (L x yz ) yz , is expected to be symmetric and positive definite [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF]. The extra-diagonal coefficients describe how pairs of processes interact, they vanish when the two corresponding processes are assumed to be independent.

Phenomenological laws of transport processes

While the symmetry principle is a deep and general argument to motivate the uncoupling of some processes, we can use more ad-hoc assumptions to simplify the entropy production, especially the term Ψ gradient from equation ( 39) that gathers all the transport processes. Here we will assume that there is no interaction between two transport processes occurring through different parts of the system. In addition, considering a moving fluid α ∈ P, we will assume that the flow process (q α ) and the diffusive processes (J u α and (J n,i α ) i∈C ) are independent. These assumptions result in splitting the term Ψ gradient in several independently positive sub-terms:

Ψ gradient = Ψ J r + {k,l}∈I(P∪{r}) Ψ J {k,l} + α∈P Ψ J α + Ψ q α Ψ J r = -J u r • ∇ τ r ≥ 0 Ψ J {k,l} = -J u {k,l} • ∇ τ {k,l} ≥ 0 Ψ J α = -J u α • ∇ τ α - i∈C J n,i α • ∇ ν i α ≥ 0 Ψ q α = - 1 T α q α • (∇ p α -ρ α g) ≥ 0 . (40)
Under the assumptions of section 4.2, the following will establish and discuss the linear phenomenological laws modelling the fluxes involved in the mass and energy balance equations of section 3.

Fluid flow

According to the entropy production Ψ q α ≥ 0 of equation ( 40), the linear phenomenological law for the velocity of the fluid phase α ∈ P is

q α = -L q α (∇ p α -ρ α g) with L q α ≥ 0 . (41)
This equation corresponds to the Darcy-Muskat law [START_REF] Leverett | Capillary behavior in porous solids[END_REF][START_REF] Muskat | The flow of heterogeneous fluids through porous media[END_REF][START_REF] Muskat | Flow of gas-liquid mixtures through sands[END_REF] with L q α as the mobility of the fluid α. The mobility is usually split between the absolute permeability of the matrix, K, the relative permeability of the fluid in the matrix, k α , and the dynamic viscosity of the fluid, µ α :

L q α = k α µ α K .
This demonstrates that the Darcy-Muskat law is compatible with the principles of thermodynamics, as already shown by Marle in [START_REF] Marle | On macroscopic equations governing multiphase flow with diffusion and chemical reactions in porous media[END_REF].

Heat conduction in matrix

According to the entropy production Ψ J r ≥ 0 of equation ( 40), the linear phenomenological law for the heat conduction through the matrix is

J u r = -L u r ∇ τ r with L u r ≥ 0 . ( 42
)
From the definition (36) of τ r , this flux can be expressed with the gradient of temperature by

J u r = - 1 T 2 r L u r ∇ T r
which corresponds to the well-known Fourier's law with 1 T 2 r L u r as effective conductivity for the matrix.

Heat conduction in interfaces

According to the entropy production Ψ J {k,l} ≥ 0 of equation ( 40), the linear phenomenological law for the heat conduction through the interface {k, l} ∈ I(P ∪ {r}) is J u {k,l} = -L u {k,l} ∇ τ {k,l} with L u {k,l} ≥ 0 .

This contribution to the total heat flow appears in some references [START_REF] Hassanizadeh | Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries[END_REF][START_REF] Marle | On macroscopic equations governing multiphase flow with diffusion and chemical reactions in porous media[END_REF] but is generally neglected. This approximation would correspond to L u {k,l} = 0.

Matter diffusion and heat conduction in fluids

According to the entropy production Ψ J α ≥ 0 of equation ( 40), the linear phenomenological law for the diffusion and conduction through the fluid α ∈ P is

               J u α = -L uu α ∇ τ α - j∈C L un, j α ∇ ν j α J n,i α = -L nu,i α ∇ τ α - j∈C L nn,i j α ∇ ν j α ∀i ∈ C ( 43
)
where the phenomenological coefficients form a symmetric positive definite matrix:

L J α =                      L uu α • • • L un, j α • • • . . . • • • • • • • • • L nu,i α • • • L nn,i j α • • • . . . • • • • • • • • •                     
.

Like for the rock matrix, the coefficient L uu α is related to the effective conductivity of the Fourier's law by 1

T 2 α L uu
α . The diffusion of components in the fluid phase are not described by a classical Fick's law like J n,i α = D i α ∇ c i α . Instead, the equation ( 43) introduces the matrix L nn α that links all the diffusive fluxes, (J n,i α ) i∈C , to all the gradients of the chemical potentials, (∇ ν i α ) i∈C . While the Fick's law is a good approximation of some limiting cases, it does not satisfy equation [START_REF] Schweizer | Hysteresis in porous media: Modelling and analysis[END_REF] in general [START_REF] Kondepudi | Modern thermodynamics: from heat engines to dissipative structures[END_REF]. Other models, like the Maxwell-Stefan equation or some generalizations of the Fick's law [START_REF] Bird | Introductory transport phenomena[END_REF][START_REF] Taylor | Multicomponent mass transfer[END_REF], are fully compatible with the thermodynamics of multi-component mixtures and can be described by the matrix L nn α . Finally, the coefficients L nu α and L un α describe the cross-effects resulting from the interactions between heat and matter flows: the Soret effect, where heat flow drives a flow of matter; and the Dufour effect, where concentration gradients drive a heat flow.

Marle suggests in [START_REF] Marle | On macroscopic equations governing multiphase flow with diffusion and chemical reactions in porous media[END_REF] that the flux J n,i α is not simply the average over the REV of the molecular diffusion at pore scale but that it also includes dispersion effects due to the difference between convection at pore scale and its average over the REV. However the splitting assumed in equation ( 40) is not compatible with such cross-effect between q α and J n,i α . Besides, dispersion appears to be a non-linear effect that is negligible at small velocity, which makes the linear regime unable to correctly describe it. More suited assumptions are needed to properly take into account the dispersion effect.

Thermal and chemical local equilibria

According to the entropy production Ψ difference of equation [START_REF] Schlüter | Pore-scale displacement mechanisms as a source of hysteresis for two-phase flow in porous media[END_REF], and neglecting all the cross-effects, the linear phenomenological laws for local exchanges are

Λ u p,p = -r u p,p (τ p -τ p ) , r u p,p ≥ 0 ∀p, p ∈ P ∪ {r} ∪ I(P ∪ {r}) Λ n,i α,β = -r n,i α,β (ν i α -ν i β ) , r n,i α,β ≥ 0 ∀i ∈ C ∀α, β ∈ P .
Here, the phenomenological coefficients r are inversely proportional to the characteristic time of the corresponding exchange process.

When the flow through the porous media is slow enough, the exchange characteristic times are negligible (r → ∞), then the different parts of the system can be assumed to be locally at equilibrium at any time. As the phenomenological coefficients r degenerate while the differences vanish, the above linear laws become unsuitable. Instead, under the assumption of thermal and chemical local equilibria, they are replaced by the identities: τ = τ p ∀p ∈ P ∪ {r} ∪ I(P ∪ {r})

ν i = ν i α ∀α ∈ P ∀i ∈ C (44)
where τ and ν i denote the common values for the different parts of the system. By inverting the definitions [START_REF] Muskat | Flow of gas-liquid mixtures through sands[END_REF], the equilibrium temperature and the equilibrium chemical potential of component i are defined by

T = 1 -τ and µ i = ν i -τ ∀i ∈ C . ( 45
)
Since the exchange terms can no longer be retrieved, the total balance equations ( 26) and ( 32) should be preferred over the partial balance equations ( 22), ( 28), ( 29) and (30).

Interface local equilibrium

When the thermal and chemical equilibria ( 44) is assumed, thanks to the equations ( 6), ( 13) and [START_REF] Everett | A general approach to hysteresis. part 2: Development of the domain theory[END_REF], and the definitions ( 27), ( 33), ( 35) and ( 45), the time variation of the total internal energy of the porous system takes the form

∂u ∂t = T ∂s ∂t + i∈C µ i ∂n i ∂t -∆W I (46) with ∆W I ≡ α∈P p α ∂φ α ∂t - {k,l}∈I(P∪{r}) γ {k,l} ∂a {k,l} ∂t . ( 47
)
This form is similar to that of a simple homogeneous system but corrected by ∆W I . Here, p α ∂ t φ α represents the rate of mechanical work given by the fluid α ∈ P and γ {k,l} ∂ t a {k,l} the rate of mechanical work received by the interface {k, l} ∈ I(P ∪ {r}). Like so, ∆W I summarizes the exchanges of mechanical work between fluids and interfaces. Thanks to equation [START_REF] Schlüter | Pore-scale displacement mechanisms as a source of hysteresis for two-phase flow in porous media[END_REF], this term is also directly connected to the entropy production:

∆W I = T Ψ time derivative ≥ 0 .
Thus ∆W I is a dissipation term, that counts for the mechanical deformations and displacements of the interfaces. It is specific to multiphase flows since ∆W I vanishes when there is only one fluid phase (Card(P) = Card(I(P ∪ {r})) = 1) with a not deformable matrix. When ∆W I = 0, the interface transformations are reversible and the interfaces are at equilibrium throughout the entire process. When ∆W I > 0, the transformations are irreversible and the interfaces are out of equilibrium.

The assumption of local equilibrium for the interfaces means the system evolves under the constraint ∆W I = 0 at any time. Since all the transformations of the system satisfy this equality, the time derivatives of equation ( 47) can be converted into total derivatives to give the following constraint:

α∈P p α dφ α = {k,l}∈I(P∪{r}) γ {k,l} da {k,l} . (48)
Equation ( 48) is an alternative definition of the local equilibrium of the interfaces. The next section explores how it can help to understand capillary effects.

Capillarity

Capillary effects are complex, they arise from the surface interactions between, at least, three phases: the matrix and two fluids. Their thermodynamic bases have been explored in past decades [START_REF] Everett | Thermodynamics of multiphase fluids in porous media[END_REF][START_REF] Hassanizadeh | Thermodynamic basis of capillary pressure in porous media[END_REF][START_REF] Morrow | Physics and thermodynamics of capillary action in porous media[END_REF]. These works rely on the Young-Laplace equation that relates the pressure difference from either side of an interface to its curvature or, for spherical meniscus, to the contact angle. The drawback of this mechanical approach is to require some upscaling procedure to fully establish a model. Indeed, the Young-Laplace equation is essentially a pore scale relationship and the interface curvature or the contact angle have no meaning at REV scale.

The analysis developed here do not rely on such pore scale closure equation. Instead, it studies the interface dissipation ∆W I ≥ 0 defined at equation (47) to show that the capillary pressure curve (53) is a general consequence of the interface equilibrium as expressed in equation (48). This result holds for system with more than two fluids and gives constraints on what are compatible multiphase capillary curves. The case of the dissipative interface transformations is then discussed and connected to more complex behaviours like the dynamical capillary pressure or the drainage-imbibition hysteresis.

Grand potential of the fluid mixture

A useful quantity to carry out the analysis of ∆W I is the total volume grand potential of the fluid mixture,

ω f = α∈P φ α ω α + {k,l}∈P∪{r} a {k,l} ω {k,l} .
Here, the name fluid mixture designates the whole system minus the matrix part, or equivalently the subsystem composed by all the fluids and all the interfaces in the system volume. Using the equations (47), ( 12), ( 7) and ( 14), the time variation of ω f follows

∂ t ω f = -s I ∂ t T - α∈P φ α ∂ t p α -∆W I , (49)
where s I is the total volume entropy of the interfaces

s I = {k,l}∈P∪{r}
a {k,l} s {k,l} .

Like for the total internal energy in equation ( 46), ∆W I appears as a dissipative term that indicates the difference in the system evolution depending on whether the interfaces follow reversible or irreversible processes.

Interface equilibrium and capillary pressure curve

Under the assumption of local equilibrium of the interfaces, ∆W I vanishes and, thanks to equation (48), the variations of ω f follow an equation stricter than (49):

dω f = -s I dT - α∈P φ α dp α . (50)
This relation indicates that the grand potential of the fluid mixture is a function of the temperature and fluid pressures: ω f (T, (p α ) α ). The interface entropy s I and the fluid volume fractions (φ α ) α also become functions of the same variables, fully determined by the derivatives of ω f (T, (p α ) α ):

s I (T, (p α ) α ) = - ∂ω f ∂T (T, (p α ) α ) , φ α (T, (p β ) β ) = - ∂ω f ∂p α (T, (p β ) β ) ∀α ∈ P . (51)
In other words, under the interface equilibrium assumption, s I and (φ α ) α are state variables and ω f (T, (p α ) α ) is an equation of state for the fluid mixture. This equation of state is subject to some constraints. The first constraint is concavity. It is a consequence of the equation (49) applied on small irreversible transformations between two equilibrium states y 1 and y 2 , which gives

ω f (y 2 ) -ω f (y 1 ) -∂ y ω f (y 1 ) • (y 2 -y 1 ) ≈ - 2 1 ∆W I ≤ 0 .
For instance, the concavity implies that ∂φ α ∂p α ≥ 0 : increasing the pressure of one fluid against the others tends to increase its volume fraction. The equation of state ω f is also subject to the constraint

β∈P ∂ω f ∂p β (T, (p α ) α ) = -φ
where φ = α φ α is the porosity, which is constant since the matrix is non deformable. This constraint also reads as the directional derivative

f T,(p α ) α = -φ with f T,(p α ) α : x → ω f (T, (p α + x) α ) .
Which gives by integration, for any pressure value p 0 ,

ω f (T, (p α + p 0 ) α ) = ω f (T, (p α ) α ) -φp 0 . (52) 20
Using this relation with a fixed p 0 in equation (51) shows that

∀β φ β (T, (p α + p 0 ) α ) = φ β (T, (p α ) α )
which means that the volume fractions do not depend on the absolute value of the pressures but only on their relative differences.

In the two phase case, with P = {w, n}, the constraint (52) finds a simpler expression by introducing the function F : x → -ω f (T, 0, x)/φ. Equations ( 51) and (52) give the relation φF (p np w ) = φ n (T, p w , p n ) which takes a more usual form as

p n -p w = p c (S n ) with S n = φ n φ and p c = (F ) -1 . ( 53
)
This is the well known equation of the capillary pressure curve [START_REF] Brooks | Hydraulic properties of porous media[END_REF][START_REF] Leverett | Capillary behavior in porous solids[END_REF]. The functions ω f and p c appears to be equivalent in the two phase case. This provides a thermodynamic interpretation to the capillary pressure curve based on the interface equilibrium and the non deformable matrix. Moreover, the thermodynamic approach naturally extends to any number of fluids and establishes that only one function, ω f (T, (p α ) α ), is needed to capture the capillary behaviour of the multiphase system.

On dissipation in capillary processes

When the interface transformations are irreversible, the interfaces can no longer be considered at equilibrium during the system evolution. Interface equilibrium could still be reached as steady state, entropy being dissipated during the process. Or interfaces could be forced to jump between equilibrium states through irreversible transitions. In the followings, two kind of processes are considered. First, the relaxation towards equilibrium, where the system is out of, but still close to, equilibrium. In this case, the linear regime can be applied and gives a behavior of dynamic capillary pressure. Second, metastability at pore scale is discussed and linked to hysteresis behavior. In both cases, the dissipation in capillary processes invalidates the traditional capillary pressure curve (53) established under interface equilibrium and requires more elaborated relationship.

Relaxation towards equilibrium

Let ω * f , s * I and (φ * α ) α denote the value of ω f , s I and (φ α ) α at steady state. If these values appear to be functions of the state (T, (p α ) α ), then the steady state is considered to be an equilibrium state and ω * f , s * I and (φ * α ) α are state functions that follow the relation (50). In particular, the following equation holds

∂ ω * f (T, (p α ) α ) ∂t = -s * I (T, (p α ) α ) ∂T ∂t - α∈P φ * α (T, (p α ) α ) ∂p α ∂t ,
while the quantities ω f , s I and (φ α ) α satisfy equation (49). Rearranging and subtracting the two gives a new expression for the entropy production due to dissipation at interface:

-∂ t (ω f -ω * f ) -(s I -s * I )∂ t T - α∈P (φ α -φ * α )∂ t p α = ∆W I = T Ψ time derivative ≥ 0 .
The dissipation find here a form more compatible with the flow-force products discussed in section 4.2. For a two-phase flow, P = {n, w}, it becomes

-∂ t (ω f -ω * f ) -(s I -s * I )∂ t T -(φ n -φ * n )∂ t (p n -p w ) ≥ 0 .
Assuming the system to be close to equilibrium and neglecting cross-effects, it suggests the following phenomenological law

φ n = φ * n -A∂ t (p n -p w ) , A ≥ 0 .
This corresponds to an alternative form of the dynamic capillary pressure described in [START_REF] Hassanizadeh | Thermodynamic basis of capillary pressure in porous media[END_REF][START_REF] Kalaydjian | Dynamic capillary pressure curve for water/oil displacement in porous media: Theory vs. experiment[END_REF].

Metastability and hysteresis

When a system is driven by a cyclic variation of one external variable, the path it follows in the thermodynamic state space depends on the processes involved. If the transformations are reversible, the system state goes back and forth along a line of equilibrium states that connects the two extremal positions. If irreversible transformations occur, the path connecting the two extremal positions changes depending on its direction, such that, on a full cycle, the two paths form a loop. The presence of such loops is often referred as hysteresis but some distinctions are required. For instance, relaxation processes, like for the dynamic capillary pressure, cause a loop that is time-dependent: the slower is the transformation, the less it dissipates and the smaller is the amplitude of the loop. At extremely slow changes, the back and forth paths become indistinguishable: the loop disappears and the transformation can be approximated as reversible. However, there are systems that exhibit time-independent loop, like the drainage-imbibition cycle [START_REF] Morrow | Capillary equilibrium in porous materials[END_REF]. In this case, the paths are irreversible but also stable and reproducible. The irreversibility is highlighted by the existence of the scanning curves: at any point, reversing the changes on the system will not cause the system to take its previous path backward but a new path instead. All states on the paths are stable since the loop is time-independent. In particular, when stopping the external change, the system stays indefinitely in the same state.

The domain theory [START_REF] Enderby | The domain model of hysteresis. part 1.-independent domains[END_REF][START_REF] Enderby | The domain model of hysteresis. part 2.-interacting domains[END_REF][START_REF] Everett | A general approach to hysteresis. part 2: Development of the domain theory[END_REF][START_REF] Everett | A general approach to hysteresis[END_REF] proposes a general thermodynamics understanding of the latter type of hysteresis, associated with time-independent loops. Two recent overview articles [START_REF] Albers | Modeling the hysteretic behavior of the capillary pressure in partially saturated porous media: a review[END_REF][START_REF] Schweizer | Hysteresis in porous media: Modelling and analysis[END_REF], dedicated on modelling capillary hysteresis in porous media, develop its historical context and its extensions. The theory is general in the sense that it is not restricted to porous media: the hysteresis phenomenon is explained as the result of a large assembly of small metastable switches. Here, a system is said to be metastable for some given control variables when it possesses two, or more, stable equilibrium states at the same values of control variables. This is only possible if the system have internal, uncontrolled, variables. When the control variables change, the equilibrium states change too, each following its own branch. The metastable system can evolve reversibly on an equilibrium branch as long as the equilibrium is stable. If the equilibrium of a branch become unstable, the branch stops. A switch is a metastable system whose at least one equilibrium branch stops. When the switch state reaches a branch stop, it rapidly jumps to the remaining stable equilibrium state through an irreversible and dissipative process. The domain theory shows that a system made of a large assembly of metastable switches, each slightly different from the others, presents the properties of hysteresis, like time-independent loops and scanning curves. The metastable switches contain elementary hysteresis, with discontinuous jumps, but the averaging over their large assembly generates the smooth, stable and reproducible loops.

For the hysteresis of the drainage-imbibition cycle, several mechanisms at pore scale have been considered as the main sources of hysteresis. The most regarded are the "ink bottle" effect, the "rain-drop" effect and the "snap-off" effect. (i) In the "ink bottle" effect [START_REF] Miller | Physical theory for capillary flow phenomena[END_REF], the only stable states for the pores are when they are filled with a unique fluid. The fluid-fluid interfaces are blocked at the bottle-necks and when the pressure difference around an interface reaches a threshold, the pore fills up rapidly with new fluid. The pores act as metastable switches where the interfaces jump between bottle-necks. (ii) The "rain drop" effect [START_REF] De Gennes | Wetting: statics and dynamics[END_REF] designates the contact angle hysteresis that appears at the solid-fluid-fluid triple line. It is named after the shape taken by a drop flowing on a plate: the contact angles are different at the front and at the rear. Here, the smooth and stable hysteresis loop is already present at pore scale, thus it is expected to be also observed at large scale. The sources of the contact angle hysteresis are found in metastability in the triple line, caused for instance by surface roughness or adsorption. (iii) The "snap-off" effect [START_REF] Schlüter | Pore-scale displacement mechanisms as a source of hysteresis for two-phase flow in porous media[END_REF][START_REF] Singh | Dynamics of snap-off and pore-filling events during two-phase fluid flow in permeable media[END_REF][START_REF] Wardlaw | The effects of geometry, wettability, viscosity and interfacial tension on trapping in single pore-throat pairs[END_REF] refers to the trapping of disconnected fluid. The effect is difficult to describe and model. For instance, the fact that fluid phases become disconnected weaken the assumption of one macroscopic pressure by phase. The metastability comes from the additional states the trapped fluid can reach.

A formulation for multiphase compositional flows

This section draws an explicit model for multiphase compositional flows in porous media from the developments above. First the existence of equations of state for the all porous system is discussed along with its consequences. Then the formulation is written as a system of partial differential equations.

The main assumptions of the model are: the matrix is inert and non deformable; no chemical reaction occurs, only change of phase of the components (vaporisation and dissolution); the porous system is always at local equilibrium; the molar and energy balance equations are given by the equations ( 26) and [START_REF] Miller | Physical theory for capillary flow phenomena[END_REF]; the fluxes are given by the equations ( 41), ( 42) and [START_REF] Taylor | Multicomponent mass transfer[END_REF]. Here, the local equilibrium of the porous system means: local equilibrium for each part of the system taken separately (see section 2.4); thermal local equilibrium (see equation ( 44)); chemical local equilibrium (see equation ( 44)); interface local equilibrium (see equation ( 48)).

Equations of state of the porous system

Considering the equations of state ω r (T ), ω f (T, (p α ) α ) and (p α (T, (µ i ) i )) α to be known, the grand potential of the porous system is the state function given by

ω(T, (µ i ) i ) = φ r ω r (T ) + ω f (T, (p α (T, (µ i ) i )) α ) .
According to equations ( 5), ( 20) and (50), its differential satisfies dω = -sdT -i∈C n i dµ i . Therefore, the total entropy of the porous is a function of the temperature and the chemical potentials: s(T, (µ i ) i ). Such an expression of the entropy is a first equation of state for the porous system, derived from the equations of state of the different parts of the system. Some other equations of state are better related to the balance equations, thanks to the change of variables ( 45):

τ = -1 T and ν i = µ i T ∀i ∈ C .
Indeed, with (τ, (ν i ) i ) the Gibbs-Duhem relation [START_REF] Bourgeat | Two-phase, partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository[END_REF] becomes

dp α = T h α dτ + T i∈C c i α dν i ∀α ∈ P and (τω) satisfies d(τω) = udτ + i∈C n i dν i .
Hence the transported quantities, (h α ) α and (c i α ) α,i , and the conserved quantities, u and (n i ) i , of the balance equations ( 26) and ( 32) become conjugate with the intensive variables (τ, (ν i ) i ). More particularly, (τω) as a function of (τ, (ν i ) i ) appears to be the Legendre transform of the total entropy of the system:

τω -τu - i∈C ν i n i = 1 T         -ω + u + i∈C µ i n i         = s with ds = -τdu - i∈C ν i dn i .
This provides a new equation of state for the system, s(u, (n i ) i ), expressed in term of total averaged variables rather than intensive variables like with s(T, (µ i ) i )) or (τω)(τ, (ν i ) i ). An other consequence is the existence of a bijection between the sets of variables (τ, (ν i ) i ) and (u, (n i ) i ) when the equations of state are smooth. While the equations of state must be convex, they are not necessarily smooth. In particular, ω f will have non differentiable points when the associated capillary pressure curve accepts a range of saturation values for the same capillary pressure. In this case, the total derivatives should be replaced by subdifferentials, and the bijection between (τ, (ν i ) i ) and (u, (n i ) i ) should be replaced by the fact that (τ, (ν i ) i , u, (n i ) i ) belong to a (1 + Card(C))-dimensional manifold.

Elliptic-parabolic differential equations

Let n = 1+Card(C) denote the number of unknowns of the problem and d = 1, 2 or 3 the space dimension. Now let v ∈ R n be the vector of primary variables defined by v = (v i ) i∈{0}∪C with

v i =        τ if i = 0 ν i if i ∈ C .
The balance equations for the internal energy of the system and the total number of moles of each component can be written as a system of n elliptic-parabolic partial differential equations:

∂B i (v) ∂t -div         j∈{0}∪C A i, j (v) ∇ v j -C i (v)         = F i (v, ∇ v) ∀i ∈ {0} ∪ C . (54)
Here B(v) ∈ R n and F(v, ∇ v) ∈ R n are vectors of scalars, C(v) ∈ R d n is a vector of vectors and A(v) ∈ R d×d n×n is a matrix of tensor. These coefficients are defined in the following, where the dependency of the physical quantities to the primary variables, v, is omitted to lighten the reading.

B(v) is given by the derivatives of the convex function (τω) with respect to v:

B i (v) = ∂(τω) ∂v i ∀i ∈ {0} ∪ C . (55)
F(v, ∇ v) is given by

F i (v, ∇ v) =        -g • α∈P ρ α L q α i∈{0}∪C ∂p α ∂v i ∇ v i -ρ α g if i = 0 0 if i ∈ C . (56)
C(v) is given by

C i (v) = α∈P 1 T ∂p α ∂v i ρ α L q α g ∀i ∈ {0} ∪ C . (57)
A(v) is given by A i, j (v) = A u i, j (v) + α∈P L J α,i, j + 1 T ∂p α ∂v i ∂p α ∂v j L q α ∀i, j ∈ {0} ∪ C (58)

where the A u i, j (v) are tensors given by

A u i, j (v) =        L u r if i = j = 0 0 if i ∈ C or j ∈ C .
The formulation (54)-( 55)-( 56)-( 57)-( 58) is fully characterized by the following inputs: the equations of state of the matrix, ω r (T ), of the fluids, p α (T, (µ i ) i ) with α ∈ P, and of the fluid mixture, ω f (T, (p α ) α ) ; and the phenomenological coefficients for the conduction in the matrix, L u r (v), for the conduction-diffusion in the fluids, L J α (v), and for the fluid velocities, L q α (v). All other quantities involved in the formulation can be retrieved from these inputs and the values of the primary variables v.

The formulation uses a persistent set of primary variables that does not depend on which phases are locally present. This feature has already been obtained for twophase flows [START_REF] Abadpour | Method of negative saturations for modeling two-phase compositional flow with oversaturated zones[END_REF][START_REF] Bourgeat | Two-phase, partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository[END_REF][START_REF] Bourgeat | On persistent primary variables for numerical modeling of gas migration in a nuclear waste repository[END_REF][START_REF] Jaffré | Henry'law and gas phase disappearance[END_REF][START_REF] Neumann | Modeling and simulation of twophase two-component flow with disappearing nonwetting phase[END_REF], but here it holds for an arbitrary numbers of phases and components, and no preferential phase is assumed. Besides, the system of equations ( 54) is well structured since the accumulation term B(v) derives from a convex potential, (τω)(v), and the elliptic operator A(v) is symmetric and positive. Furthermore, the entropy equation ( 34) is directly retrieved by summing the products of the equations (54) with the (v i ) i . Among the others formulations [START_REF] Abadpour | Method of negative saturations for modeling two-phase compositional flow with oversaturated zones[END_REF][START_REF] Bourgeat | Two-phase, partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository[END_REF][START_REF] Bourgeat | On persistent primary variables for numerical modeling of gas migration in a nuclear waste repository[END_REF][START_REF] Class | Numerical simulation of non-isothermal multiphase multicomponent processes in porous media.: 1. an efficient solution technique[END_REF][START_REF] Coats | An equation of state compositional model[END_REF][START_REF] Forsyth | A two-phase, two-component model for natural convection in a porous medium[END_REF][START_REF] Jaffré | Henry'law and gas phase disappearance[END_REF][START_REF] Neumann | Modeling and simulation of twophase two-component flow with disappearing nonwetting phase[END_REF] that describe a comparable range of physical behaviors, the majority lacks this level of generality and none of them gather such structuring properties.

Concluding remarks

This paper applies the principles of the thermodynamics of irreversible processes to the representative elementary volume of a porous medium filled with multiple compositional fluids in order to model heat and mass transfers. A special attention has been taken to treat the interfaces in a consistent way. The physical assumptions made are usual and the resulting model does not present new macroscopic behaviors. Nevertheless, the final formulation exhibits less common, but useful, properties such as a persistent set of primary variables, a symmetric structure and an associated entropy equation.

One benefit of the thermodynamic approach is to avoid assuming specific processes at the microscopic scale, only that the system is close to equilibrium. This general framework turns out to be sufficient to justify macroscopic behaviors like the Darcy-Muskat law and the capillary pressure curve. The proposed explanation of capillarity in particular contrasts with the existing works that all rely on the Young-Laplace equation at pore scale. Moreover, it provides new and strong constraints on capillary pressure curve for three or more phases.

The physical assumptions made in this paper should be relaxed to include additional processes. For instance, taking the reactive chemistry into account could be of peculiar interest because chemistry has already a strong thermodynamic basis and the present work gives a thermodynamic understanding of flows in porous media. Formulating reactive transport on this common ground could help in properly solving this kind of coupling.

  indicates that the internal energy of the interface {k, l} is a function of the form U {k,l} (S {k,l} , A {k,l} ) and that its derivatives are the temperature T {k,l} and the surface tension γ {k,l} of the interface {k, l}. According to the assumptions, all the interfaces between phases k and l in the system volume V have the same temperature and the same surface tension. Applying the Euler's homogeneous function theorem gives U {k,l} = T {k,l} S {k,l} + γ {k,l} A {k,l} ∀ {k, l} ∈ I(P ∪ {r})

  a {k,l} u {k,l} . According to equation[START_REF] Cheng | Heat transfer in geothermal systems[END_REF], we have a {k,l} u {k,l} = a {k,l} T {k,l} s {k,l} + a {k,l} γ {k,l} ∀ {k, l} ∈ I(P ∪ {r}) , the differential of this relation gives d(a {k,l} u {k,l} ) = T {k,l} d(a {k,l} s {k,l} ) + γ {k,l} da {k,l} + a {k,l} s {k,l} dT {k,l} + a {k,l} dγ {k,l} , Similarly, the amount of grand potential of the interface {k, l} per unit of system volume is Ω {k,l} /V = a {k,l} ω {k,l} . According to equation[START_REF] Coats | An equation of state compositional model[END_REF], its differential satisfies d(a {k,l} ω {k,l} ) = d(a {k,l} u {k,l} )d(a {k,l} T {k,l} s {k,l} ) ,

	which, using equation (12), reduces to
	d(a {k,l} u {k,l} ) = T {k,l} d(a {k,l} s {k,l} ) + γ {k,l} da {k,l} ∀ {k, l} ∈ I(P ∪ {r}) .	(13)
	which, using equation (13), reduces to
	d(a

{k,l} ω {k,l} ) = -a {k,l} s {k,l} dT {k,l} + γ {k,l} da {k,l} ∀ {k, l} ∈ I(P ∪ {r}) .