Joubine Aghili

Jean-Raynald De Dreuzy

Laurent Trenty

Roland Masson

A hybrid

Keywords: Discrete Fracture Models, two-phase Darcy flow, Fickian diffusion, phase transitions, finite volume scheme, two-point flux approximation

Introduction

Flow and transport in fractured porous media are key processes in many subsurface applications such as geothermal energy, CO 2 sequestration, hydrogeology, or geological storage. In this work, we consider the class of model representing the fractures as co-dimension one surfaces coupled with the surrounding matrix. These so called hybrid-dimensional or Discrete Fracture Matrix (DFM) models are obtained by averaging both the unknowns and the equations in the fracture width and by imposing appropriate transmission conditions at both sides of the matrix fracture (mf) interfaces. Given the high geometrical complexity of real-life fracture networks, the main advantages of these hybrid-dimensional models compared with full-dimensional ones are to facilitate the mesh generation and the discretization of the model, and to reduce the computational cost of the resulting schemes. This type of reduced models has been the object of intensive researches over the last 15 years due to the ubiquity of fractures in geology and their considerable impact on the flow, transport of mass and energy in porous media.

For single phase Darcy flow, two types of hybrid-dimensional models are classically distinguished. First, in the case of fractures acting as drains, the continuity of the pressure can be assumed at mf interfaces leading to the so-called continuous pressure models. Continuous pressure models have been first introduced in [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF] and later analyzed for general network of planar fractures in [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media[END_REF]. To account for fractures acting both as drains or barriers, one must allow the discontinuity of the pressure at mf interfaces leading to discontinuous pressure models with additional unknowns on both sides of the mf interfaces. Discontinuous pressure models have been introduced in [START_REF] Granet | A single-phase flow simulation of fractured reservoir using a discrete representation of fractures[END_REF][START_REF] Flauraud | Domain decomposition for an asymptotic geological fault modeling[END_REF][START_REF] Karimi-Fard | An efficient discrete-fracture model applicable for general-purpose reservoir simulators[END_REF], then generalized as a one parameter class of models in [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF] and investigated by many authors with different types of discretizations [START_REF] Angot | Asymptotic and numerical modelling of flows in fractured porous media[END_REF][START_REF] Tunc | A model for conductive faults with non-matching grids[END_REF][START_REF] Sandve | An efficient multi-point flux approximation method for Discrete Fracture-Matrix simulations[END_REF][START_REF] Fumagalli | A reduced model for flow and transport in fractured porous media with non-matching grids[END_REF][START_REF] Schwenck | Dimensionally reduced flow models in fractured porous media: crossings and boundaries[END_REF][START_REF] Ahmed | Threedimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF][START_REF] Brenner | Gradient discretization of hybriddimensional Darcy flow in fractured porous media with discontinuous pressures at matrix-fracture interfaces[END_REF][START_REF] Antonietti | Mimetic finite difference approximation of flows in fractured porous media[END_REF][START_REF] Nordbotten | Unified approach to discretization of flow in fractured porous media[END_REF]. Let us also refer to [START_REF] Hennicker | Modeling and Analysis of the Coupling in Discrete Fracture Matrix models[END_REF] for a derivation and error analysis of these reduced models using a Fourier approach.

In the case of two-phase Darcy flow, a natural extension of continuous pressure models is considered in [START_REF] Bogdanov | Two-phase flow through fractured porous media[END_REF][START_REF] Reichenberger | A mixed-dimensional finite volume method for two-phase flow in fractured porous media[END_REF][START_REF] Monteagudo | Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects[END_REF][START_REF] Matthai | Finite element -node-centered finite-volume two-phase-flow experiments with fractured rock represented by unstructured hybrid-element meshes[END_REF][START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF][START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal darcy flows in fractured porous media[END_REF][START_REF] Brenner | Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF][START_REF] Brenner | Vertex Approximate Gradient Discretization preserving positivity for two-phase Darcy flows in heterogeneous porous media[END_REF] for fractures acting as drains. It is obtained by assuming that both phase pressures are continuous at mf interfaces. However, as noticed e.g. in [START_REF] Brenner | Hybrid dimensional modelling of two-phase flow through fractured with enhanced matrix fracture transmission conditions[END_REF][START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF], this assumption can lead to inaccurate approximations of the mf fluxes due to the inconsistency of matching the mf interfaces capillary pressure with the one inside the fracture. An alternative hybrid-dimensional model is widely used in the case of cell centered discretizations such as in [START_REF] Granet | A two-phase flow simulation of a fractured reservoir using a new fissure element method[END_REF][START_REF] Karimi-Fard | An efficient discrete-fracture model applicable for general-purpose reservoir simulators[END_REF] for Two-Point Flux Approximations (TPFA) and in [START_REF] Gläser | A discrete fracture model for twophase flow in fractured porous media[END_REF][START_REF] Gläser | A hybrid-dimensional discrete fracture model for non-isothermal two-phase flow in fractured porous media[END_REF] for Multi-point Flux Approximations (MPFA). Let us also refer to [START_REF] Brenner | Hybrid dimensional modelling of two-phase flow through fractured with enhanced matrix fracture transmission conditions[END_REF] for an extension to the case of the Vertex Approximate Gradient (VAG) scheme. In these approaches, the transmissivities are first computed for a linear single phase Darcy flow model with elimination of the mf interface (or sub-interface for MPFA) pressures from the flux continuity equations. Then, the mf fluxes are obtained combining these single phase Darcy fluxes with an upwind approximation of the mobilities jumping over the mf interfaces i.e. using cell and fracture face saturations and not mf interface saturations. Using the terminology introduced in [START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF], this type of model is termed mf linear f-upwind hybrid-dimensional model. As investigated in [START_REF] Brenner | Hybrid dimensional modelling of two-phase flow through fractured with enhanced matrix fracture transmission conditions[END_REF][START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF], its main drawback is that it does not capture accurately the saturation jumps at mf interfaces, leading typically to an over-estimated leakage for fracture acting as drains, especially on coarse meshes. Let us also refer to [START_REF] Hoteit | An efficient numerical model for incompressible twophase flow in fractured media[END_REF][START_REF] Moortgat | Higher-order compositional modeling of three-phase flow in 3d fractured porous media based on cross-flow equilibrium[END_REF][START_REF] Zidane | Fracture-cross-flow equilibrium in compositional twophase reservoir simulation[END_REF] for related approaches based on the cross flow equilibrium concept and to [START_REF] Jiang | An improved projection-based embedded discrete fracture model (pedfm) for multiphase flow in fractured reservoirs[END_REF] -and references there-in -for non-conforming approximations using projection based Embedded Discrete Fracture models. To obtain a more accurate approximation of the mf fluxes one must turn to hybrid-dimensional models keeping both the pressures and the saturations interfacial unknowns. Then, these models write the half two-phase Darcy flux continuity equations on both sides of the mf interface using a TPFA of the half fluxes on the fracture side. It leads to the so-called mf nonlinear hybrid-dimensional models [START_REF] Jaffré | A discrete fracture model for two-phase flow with matrix-fracture interaction[END_REF][START_REF] Droniou | Numerical analysis of a two-phase flow discrete fracture model[END_REF][START_REF] Brenner | Hybrid dimensional modelling of two-phase flow through fractured with enhanced matrix fracture transmission conditions[END_REF][START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF][START_REF] Alali | Finite-Volume Simulation of Capillary-Dominated Flow in Matrix-Fracture Systems using Interface Conditions[END_REF][START_REF] Brenner | Total-velocity-based finite volume discretization of two-phase Darcy flow in highly heterogeneous media with discontinuous capillary pressure[END_REF] according to the terminology used in [START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF].

The main objective of this paper is to extend these types of mf nonlinear hybriddimensional models to compositional two-phase Darcy flow accounting for phase transitions and Fickian diffusion. The transmission conditions are designed to be consistent with the physical processes at mf interfaces. They account in particular for the saturation jump induced by the different rock types, for the Fickian diffusion in the fracture width, as well as for the thermodynamical equilibrium. They are based on flux continuity equations for each component using a TPFA of the fluxes in the fracture width. The saturation jumps at mf interfaces are captured for general capillary pressure curves thanks to a parametrization of the matrix and fracture capillary pressure graphs as introduced in [START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF][START_REF] Brenner | Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF]. The thermodynamical equilibrium is formulated at mf interfaces using complementary constraints and taking into account the saturation jumps.

Our choice of the formulation of the compositional model is based on [START_REF] Beaude | Non-isothermal compositional liquid gas darcy flow: formulation, soil-atmosphere boundary condition and application to high-energy geothermal simulations[END_REF]. Like the Coats' formulation [START_REF] Coats | Implicit compositional simulation of single-porosity and dual-porosity reservoirs[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] Eymard | Vertex-centred discretization of multiphase compositional Darcy flows on general meshes[END_REF][START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal darcy flows in fractured porous media[END_REF], it uses the phase pressures, phase saturations and phase molar fractions as set of principal unknowns. This is a convenient choice since all the physical laws can be directly expressed using subsets of this set of variables. In order to avoid the switch of variables, this choice of the principal unknowns is combined with an extension of the phase molar fractions of an absent phase by the molar fractions at thermodynamic equilibrium with the present phase. It results that the set of principal unknowns does not depend on the set of present phases leading to a so-called persistent variable formulation (see [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] Marchand | Fully coupled generalized hybrid-mixed finite element approximation of two-phase two-component flow in porous media. part i: formulation and properties of the mathematical model[END_REF][START_REF] Angelini | Finite volume approximation of a diffusion-dissolution model and application to nuclear waste storage[END_REF][START_REF] Masson | Formulations of two phase liquid gas compositional Darcy flows with phase transitions[END_REF] for other examples of this class of formulations). Moreover, the phase transitions can be expressed as complementary constraints as in [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF] which means that the nonlinear systems can be solved using semi-smoothed Newton techniques such as the Newton-min algorithm [START_REF] Kräutle | The semi-smooth newton method for multicomponent reactive transport with minerals[END_REF][START_REF] Gharbia | Gas phase appearance and disappearance as a problem with complementarity constraints[END_REF]. The extension to DFM models is based on pressure, saturation and molar fractions unknowns both in the matrix and in the fracture network.

The discretization of the hybrid-dimensional model is based on a TPFA combined with a phase potential upwinding of the mobilities and molar fractions. To improve the robustness of the nonlinear solver, porous volumes are added at mf interfaces by distribution of the neighboring cell and fracture face porous volumes to the mf interfaces. The same idea is applied at edges shared by at least three fractures at which a porous volume is defined by distribution of the neighboring fracture face porous volumes to that edge. To further improve the robustness, a regularization of the harmonic mean of the saturations in the definition of the Fickian diffusion fluxes is applied.

The DFM compositional model is validated by comparison with a reference equi-dimensional model using Cartesian geometries in order to be able to refine the mesh in the fracture width, even for very small apertures. Two types of test cases are considered both using a fluid system with liquid and gas phases defined as mixtures of air and water components. The first test case simulates the gas injection in an initially pure water-liquid saturated DFM with either small or high dissolution of the air component in the liquid phase. The second test case models the desaturation by suction at the interface between a low permeable and fractured porous medium and the atmosphere. The effect of the Fickian diffusion on the gas transfer at mf interfaces is investigated. To enhance the interest of the mf nonlinear model, it is also compared to the usual mf linear f-upwind model which does not incorporate mf interface saturation unknowns.

Finally, the mf nonlinear hybrid-dimensional compositional two-phase Darcy flow model is applied on a challenging test case modelling the desaturation by suction of a low permeable rock with open fractures. The set up of this test case is based on Andra nuclear waste storage prototype facility in a Callovo-Oxfordian argilite. The DFM is based on a simplified 2D version of the fracture network taking place at the interface between the exacavated tunnels and the storage rock.

The remaining of this article is outlined as follows. Section 2 describes the mf nonlinear hybrid-dimensional model using the formulation of the compositional two-phase Darcy flow introduced in [START_REF] Beaude | Non-isothermal compositional liquid gas darcy flow: formulation, soil-atmosphere boundary condition and application to high-energy geothermal simulations[END_REF]. Section 3 details the TPFA discretization of the hybrid-dimensional model including the definition of the porous volumes, the accumulation terms, the Darcy and Fick fluxes, as well as the thermodynamical equilibrium closure equations based on the complementary constraints. Then, Section 4 gathers numerical experiments including the comparisons with the reference equi-dimensional model and with the mf linear f-upwind hybrid-dimensional model, as well as the application test case.

The continuous model

We consider a compositional two-phase flow in a Discrete Fracture Matrix (DFM) model accounting for the flow and transport in the matrix coupled to the flow and transport in the network of fractures represented as co-dimension one planar surfaces.

Let Ω be an open bounded domain of R d , d = 2, 3 assumed to be polyhedral for d = 3 (and polygonal for d = 2). The spatial dimension will be fixed to d = 3 when it needs to be specified, for instance in the naming of the geometrical objects or for the space discretization in the next section. The adaptations to the case d = 2 are straightforward. Let Γ = i∈I Γ i denotes the network of fractures Γ i ⊂ Ω, i ∈ I, such that each Γ i is a planar polygonal simply connected open domain included in some plane of R d . Without restriction of generality, we will assume that the fractures may intersect exclusively at their boundaries (see Figure 1), that is for any i, j ∈ I, i = j one has

Γ i ∩ Γ j = ∅, but not necessarily Γ i ∩ Γ j = ∅.
The two sides of a given fracture of Γ are denoted by ± in the matrix domain, with unit normal vectors n ± oriented outward of the sides ±. We denote by γ ± the trace operators on the sides ± of Γ, by γ ± n the normal trace operators on the sides ± of Γ such that γ ± n q = q • n ± | Γ ± , and by • the normal trace jump operator on Γ defined by q = γ + n q + γ - n q.

We denote by ∇ τ the tangential gradient and by div τ the tangential divergence on the fracture network Γ.

The two-phase compositional model is defined by the set P = {nw, w} of "non-wetting" and "wetting" phases, each phase being a mixture of components i ∈ C. For each phase α ∈ P, we denote by ζ α (p α , c α ) the molar density, by ρ α (p α , c α ) the mass density, by µ α (p α , c α ) the dynamic viscosity, and by f α i (p α , c α ) the fugacity of the component i in phase α. They are functions of the phase pressure p α and of the phase component molar fractions c α = (c α i ) i∈C . Let also denote by s α the phase pore volume fraction or saturation.

In the matrix domain, φ m is the matrix porosity and Λ m the matrix permeability tensor, both possibly depending on x ∈ Ω. In the fracture network, φ f denotes the fracture porosity, Λ f is the tangential permeability tensor and d f the thickness of the fracture, all possibly depending on x ∈ Γ. The permeability tensor in the fracture network is assumed to have the fracture normal vector n + as principal direction and λ f,n denotes the corresponding fracture normal permeability.

For simplicity and since we focus on mf interfaces, it is assumed that both the matrix and fracture domains are homogeneous with respect to capillary pressure and relative permeability laws. For each phase α ∈ P and each rock type ν ∈ {m, f }, we denote by k α r,ν (s α ν ) the phase relative permeability as function of the phase saturation and by P c,ν (s nw ν ) the capillary pressure as function of the non-wetting phase saturation. Let us also denote by P c,ν (s nw ν ) the monotone graph extension of the capillary pressure graph. Following [START_REF] Beaude | Non-isothermal compositional liquid gas darcy flow: formulation, soil-atmosphere boundary condition and application to high-energy geothermal simulations[END_REF], the two-phase compositional model is formulated in terms of phase pressures, saturations and phase compositions (PSC). It is a permanent variable formulation based on the extension of the phase component molar fractions c α for an absent phase α (such that s α = 0) by the ones at thermodynamical equilibrium with the present phase. As the widely used Coat's formulation, the PSC formulation has the advantage to be based on the physical variables. On the other hand, it avoids the cumbersome variable switch typical of the Coat's formulation with different sets of unknowns depending on the set of present phases. Moreover, thanks to the extension of the phase component molar fractions, the thermodynamical equilibrium is expressed as complementary constraints which allows the use of Newton-min nonlinear solvers to solve the resulting nonlinear systems.

The extension of the PSC formulation to DFM models, uses as primary unknowns the phase pressures p α ν , the phase saturations s α ν and the phase component molar fractions c α ν , α ∈ P, ν ∈ {m, f }.

In addition, a parametrization of the capillary pressure graphs is used at the mf interfaces to account for the different rock types and capture the saturation jump. Denoting by τ ± the new interface variable on each side ± of the fractures, with τ ± ∈ [τ 0 , τ 1 ], this parametrization introduces a saturation function S nw ν (τ ± ) for each rock type ν ∈ {m, f } as well as a capillary pressure function P c (τ ± ) which are such that

P c (τ ± ) ∈ P c,ν (S nw ν (τ ± )), ν ∈ {m, f }, dPc dτ ± (τ ± ) + ν∈{m,f } dS nw ν dτ ± (τ ± ) > 0, and
S nw ν (τ 0 ) = 0, S nw ν (τ 1 ) = 1, ν ∈ {m, f }.
The practical choice of this parametrization is based on a switch of variables technique between the matrix, fracture saturations and the capillary pressure (see the example of parametrization Figure 2 based on a switch of variable between the fracture and the matrix

a f = 2b f , a m = 5b m , τ ∈ (τ 0 , τ 1 ) with τ 0 = 0, τ 1 = 2 - 1 -e -a f bm 1 -e -am bm .
non-wetting phase saturations). This parametrization has the advantage to account for general capillary pressure curves including typically entry pressure in the matrix and vanishing capillary pressure in the fractures. It also provides a more efficient nonlinear convergence than the more usual capillary pressure or saturation variable. Let us refer to [START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF][START_REF] Brenner | Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF] for details. In the following we set

S w ν (τ ± ) = 1 -S nw ν (τ ± ).
The set of equations couples the molar conservation of each component in the matrix domain and in the fracture network averaged over the width of the fracture

φ m ∂ t n i,m + divq i,m = 0, i ∈ C, φ f d f ∂ t n i,f + div τ (d f q i,f ) -q i,m = 0, i ∈ C, (1)
to the capillary relation and the sum to one of the saturations

p c,ν = p nw ν -p w ν = P c,ν (s nw ν ), s nw ν + s w ν = 1, (2)
and to the thermodynamical equilibrium

       f nw i (p nw ν , c nw ν ) = f w i (p w ν , c w ν ), i ∈ C, min 1 - i∈C c α i,ν , s α ν = 0, α ∈ P, (3)
for ν ∈ {m, f }. In the above conservation equations, the component number of moles per unit volume is defined by

n i,ν = α∈P ζ α (p α ν , c α ν )s α ν c α i,ν (4)
and the component molar flux by

q i,ν = α∈P (c α i,ν V α ν + F α i,ν ). (5)
Using the convenient notation ∇ m = ∇ and ∇ f = ∇ τ , the phase Darcy flux V α ν is obtained from the generalized Darcy law such that

V α ν = - ζ α (p α ν , c α ν )k α r,ν (s α ν ) µ α (p α ν , c α ν ) Λ ν ∇ ν (p α ν -ρ α (p α ν , c α ν )g ν ) ( 6
)
where g is the gravity acceleration vector and g m = g, g f = g -(g • n + )n + . The Fickian diffusion flux F α i,ν is defined by

F α i,ν = -φ ν s α ν ζ α (p α ν , c α ν ) D α T 2 ν ∇ ν c α i,ν , (7)
where D α is the Fickian diffusion coefficient in phase α assumed to be the same for all components, and T ν is the tortuosity of medium ν ∈ {m, f }.

Figure 3: (Left): example of a 2D DFM model with the normal vectors n ± at both sides of a fracture, the matrix phase pressure, saturation and molar fractions p α m , s α m , c α m , the fracture phase pressure, saturation and molar fractions p α f , s α f , c α f , the matrix component flux q i,m and the fracture tangential component flux q i,f . (Right): zoom at the interface with artificially enlarged fracture illustrating the transmission condition γ + n q i,m + q + i,f = 0 on the + side of the interface Γ with γ + n the normal trace operator on the + side of Γ, q + i,f

the component normal flux on the + side of the fracture, and τ + the parameter variable on the + side. Note also that the interface saturation S α m (τ + ) on the matrix side differs from the interface saturation S α f (τ + ) on the fracture side due to different rock types.

The system of equations ( 1)-( 2)-( 3) is complemented by transmission conditions at mf interfaces. On both sides of the interface, they account for the component normal flux conservations

     γ ± n q i,m + q ± i,f = 0, q ± i,f = α∈P (V α,± i,f + F α,± i,f ), i ∈ C. (8)
combined with the continuity of the matrix pressure up to the mf interface for phases flowing on both sides of the interface and with

γ ± s nw m = S nw m (τ ± ), (9)
relating the mf interface variable τ ± to the trace of the non-wetting saturation on the matrix side. The normal Darcy and Fickian diffusion fluxes on the fracture sides are defined by a two-point approximation between the fracture and each interface ±. To be more specific, the approximation

F α,± i,f
of the normal Fickian diffusion flux is based on a twopoint approximation of the normal derivative combined with an arithmetic averaging of the molar density and an harmonic averaging of the saturation. This harmonic averaging ensures that the Fickian diffusion vanishes if the phase disappear either in the fracture or at the interface on the fracture side

F α,± i,f = ζ α (p α f , c α f ) + ζ α (γ ± p α m , γ ± c α m ) 2 s α f S α f (τ ± ) s α f + S α f (τ ± ) φ f D α T 2 f (c α i,f -γ ± c α i,m ) d f /2 . ( 10
)
The approximation of the Darcy flux is based on the two-point approximation of the normal derivative with an arithmetic averaging of the mass density combined with a phase potential upwind approximation of the mobilities between the fracture and the interface ± on the fracture side

V α,± f = λ f,n (p α f -γ ± p α m ) d f /2 + ρ α (p α f , c α f ) + ρ α (γ ± p α m , γ ± c α m ) 2 g f • n ± , (11)
V α,± i,f = c α i,f ζ α µ α (p α f , c α f ) k α r,f (s α f ) (V α,± f ) + -γ ± c α i,m ζ α µ α (γ ± p α m , γ ± c α m ) k α r,f (S α f (τ ± )) (V α,± f ) -, (12)
where we have used the notations (a) + = max(0, a) and (a) -= -(-a) + for all a ∈ R. Note that the saturation at the interface Γ ± is discontinuous with value S α f (τ ± ) on the fracture side and S α m (τ ± ) on the matrix side (see Fig. 3). It results that the thermodynamical equilibrium (3) must be adapted at the interface such that

                 f nw i (γ ± p nw m , γ ± c nw m ) = f w i (γ ± p w m , γ ± c w m ), i ∈ C, min 1 - i∈C γ ± c nw i,m , τ ± -τ 0 = 0, min 1 - i∈C γ ± c w i,m , τ 1 -τ ± = 0, (13)
meaning that the phase disappears only if the largest phase saturation between those on the matrix and fracture sides vanishes.

Finally, normal flux conservation is also imposed at fracture intersections and we assume for simplicity that the phase component molar fractions and the phase pressures are continuous at fracture intersections.

3 The discrete model

Space discretization

We will denote by M the set of disjoint open polyhedral cells, by F the set of faces and by E the set of edges of the mesh. For each cell K ∈ M, we denote by F K ⊂ F the set of its faces and for each face σ ∈ F, we denote by E σ the set of its edges. In addition, the set of cells sharing a given face σ ∈ F is denoted by M σ , and the set of faces sharing a given edge e ∈ E is denoted by F e . The mesh is supposed to be conforming w.r.t. the fracture network Γ in the sense that, for each i ∈ I, there exists a subset

F Γi of F such that Γ i = σ∈FΓ i σ.
Let us denote by F Γ the subset i∈I F Γi of fracture faces and let us set

E Γ = σ∈FΓ E σ .
For all K ∈ M and σ ∈ F K , we denote by n Kσ the unit outward normal vector to σ w.r.t K, and similarly, for all σ ∈ F Γ and e ∈ E Γ , we denote by n σe the unit outward normal vector to e w.r.t σ. Finally, we denote by |K|, K ∈ M, the d dimensional measure of K, by |σ|, σ ∈ F the d -1 dimensional measure of σ and by |e|, e ∈ E Γ the d -2 dimensional measure of e.

The degrees of freedom of the discretization scheme are associated with the "cell centres" x K ∈ K, the "centres" of the fracture faces x σ ∈ σ, and the "centres" x e ∈ e of certain edges of the mesh. More precisely, in order to deal with the intersection of multiple fractures at a given edge e ∈ E Γ we associate to it a degree of freedom and a certain porous volume. We will denote by E the subset of E Γ such that for all e ∈ E the set F Γ ∩ F e has at least 3 elements. The set of control volumes also include the mf interfaces defined by the set

I = {(K, σ) | σ ∈ F Γ , K ∈ M σ }.
For brevity, (K, σ) will be simply denoted Kσ. The set of degrees of freedom is defined formally by

D = {K ∈ M} ∪ {σ ∈ F Γ } ∪ {Kσ ∈ I} ∪ {e ∈ E } ,
see Fig. 4 for an illustration on a 2D triangular mesh. The consistency of the two-point flux approximation is based on the orthogonality assumption in the sense that x K x σ ⊥ σ for K ∈ M σ , σ ∈ F K , and x σ x e ⊥ e for σ ∈ F Γ , e ∈ E σ . It is also assumed that the permeability tensors Λ ν are isotropic and we set

Λ K I = 1 |K| K Λ m (x)dx, Λ σ I = 1 |σ| σ Λ f (x)dσ(x), λ f,n,σ = 1 |σ| σ λ f,n (x)dσ(x).
The anisotropy of the permeability tensors can also be dealt with if the directions of anisotropy are aligned with the face normal vectors n Kσ for Λ K and with the edge normal vectors n σe in the face plane for Λ σ .

Numerical scheme

Let t 0 = 0 < t 1 < • • • < t N = T a discretization of the time interval [0, T ] with T > 0 the final time, and let us set ∆t n = t n -t n-1 for n = 1, • • • , N . The superscript n will be used in the following to denote the value at time t n . The discretization of equations ( 1)-( 12) is based on a fully implicit time integration in order to avoid strong time step restrictions due to the typically high velocities in the fractures and to cope with the strong nonlinear coupling at mf interfaces. It is combined with a finitevolume method in space using two-point flux approximation. The finite volume scheme incorporates porous volumes at mf interfaces and at fracture intersections e ∈ E in order to avoid possible degeneracies of the nonlinear and linear systems and to improve the nonlinear convergence. The discrete physical unknowns are matching the ones from the PSC model, namely the pressures {p α κ } κ∈D , the saturations {s α κ } κ∈D\I , the phase component molar fractions {c α κ } κ∈D , and the mf interface variables {τ κ } κ∈I , at each time step n ∈ {1, . . . , N }. The closure equations (2) are used to eliminate wetting phase pressure and saturation choosing the non-wetting phase pressure p nw κ , κ ∈ D as reference pressure and either the non-wetting phase saturation s nw κ for κ ∈ D \ I or the interface variable τ κ for κ ∈ I. It results that we implicitly set in the following

               p w κ = p nw κ -P c,m (s nw κ ) for κ ∈ M, p w κ = p nw κ -P c,f (s nw κ ) for κ ∈ F ∪ E , p w κ = p nw κ -P c (τ κ ) for κ ∈ I, s w κ = 1 -s nw κ for κ ∈ D \ I, 1 -s w K,Kσ = s nw K,Kσ = S nw m (τ Kσ ) for Kσ ∈ I, 1 -s w σ,Kσ = s nw σ,Kσ = S nw f (τ Kσ )
for Kσ ∈ I.

The set of physical primary unknowns at each control volume κ is denoted by X κ with

X κ = (p nw κ , s nw κ , c α κ , α ∈ P) for κ ∈ D \ I, (p nw κ , τ κ , c α κ , α ∈ P) for κ ∈ I,
and X D = {X κ , κ ∈ D} is the full set of primary unknowns. Let us denote by n i,κ the number of moles of component i in the control volume κ ∈ D, and by q i,κ,κ the discrete component normal flux connecting the control volumes κ and κ , oriented outward to κ. The definitions of these number of moles and fluxes are detailed in the next subsections 3.2.1 and 3.2.2. The system of discrete equations couples at each time step n the discrete molar conservation equations defined for i ∈ C by • for all cells K ∈ M

n n i,K -n n-1 i,K ∆t n + σ∈F K \FΓ q n i,K,σ + σ∈F K ∩FΓ q n i,K,Kσ = 0, (14)
• for all fracture faces σ ∈ F Γ

n n i,σ -n n-1 i,σ ∆t n + e∈Eσ q n i,σ,e + K∈Mσ q n i,σ,Kσ = 0, (15)
• for all edges e ∈ E

n n i,e -n n-1 i,e ∆t n - σ∈Fe∩FΓ q n i,σ,e = 0, (16)
• for all interfaces Kσ ∈ I

n n i,Kσ -n n-1 i,Kσ ∆t n -q n i,K,Kσ -q n i,σ,Kσ = 0, (17)
and the thermodynamical equilibrium such that

                             f nw i (p nw,n κ , c nw,n κ ) = f w i (p w,n κ , c w,n κ ), i ∈ C, κ ∈ D, min 1 - i∈C c α,n i,κ , s α,n κ = 0, α ∈ P, κ ∈ D \ I, min 1 - i∈C c nw,n i,κ , τ n κ -τ 0 = 0, κ ∈ I, min 1 - i∈C c w,n i,κ , τ 1 -τ n κ = 0, κ ∈ I. (18)
It is easy to check that the number of physical unknowns per control volume matches with the number of equations per control volume and is equal to 2 + 2#C. A Newton-min algorithm is used to solve the nonlinear system at each time step of the simulation. The Jacobian is reduced without additional fill-in to #C unknowns and equations per control volume by elimination of the #C + 2 thermodynamical equilibrium equations [START_REF] Reichenberger | A mixed-dimensional finite volume method for two-phase flow in fractured porous media[END_REF] and of #C + 2 unknowns. The choice of these eliminated unknowns depends on the active constraints matching with the set of present phases in our implementation of the Newton-min algorithm (see [START_REF] Beaude | Non-isothermal compositional liquid gas darcy flow: formulation, soil-atmosphere boundary condition and application to high-energy geothermal simulations[END_REF] for details).

The Jacobian is further reduced by Schur elimination of the interface unknowns and of the interface conservation equations. For Kσ ∈ I, these equations connect Kσ to the cell K and to the fracture face σ such that the only additional connections in the Schur complement are between each fracture face σ and the cells K ∈ M σ . This elimination of the interface unknowns is a crucial step which considerably improves the condition number of the Jacobian.

Porous volumes and accumulations

Let us first define the porous volumes φ κ at each control volume κ ∈ D including a redistribution of neighbouring cell and fracture face (resp. neighbouring fracture faces) porous volumes to the interfaces Kσ ∈ I (resp. to the edges e ∈ E ). Let us denote by {α K,Kσ } Kσ∈I , {α σ,Kσ } Kσ∈I (resp. {α e,σ } σ∈Fe∩FΓ,e∈E ) these redistribution coefficients.

Let φK = K φ m (x)dx and φσ = σ φ f (x)d f (x)dσ(x), then we set φ K = φK - σ∈F K ∩FΓ α K,Kσ φK , K ∈ M, φ σ = φσ - K∈Mσ α σ,Kσ φσ - e∈E ∩Eσ α e,σ φσ , σ ∈ F Γ , φ K,Kσ = α K,Kσ φK , Kσ ∈ I, φ σ,Kσ = α σ,Kσ φσ , Kσ ∈ I, φ e = σ∈Fe∩FΓ α e,σ φσ , e ∈ E . (19)
The coefficients α K,Kσ and α σ,Kσ must be strictly positive to avoid the singularity of the Jacobian matrix which could occur otherwise when the phase permeability vanishes. If they are chosen small enough their impact on the solution is negligible. In practice α σ,Kσ = 0.1 is a good choice. If the matrix acts as a barrier α K,Kσ = 0.1 is also fine thanks to the saturation jump. For a matrix acting as a drain α K,Kσ = 0.01 is better.

On the other hand, the coefficients α e,σ are chosen roughly speaking to equidistribute the porous volume at the edge with those at the neighbouring fracture faces.

Given the porous volumes, the number of moles per control volume are defined by

n i,κ = φ κ α∈P ζ α (p α κ , c α κ )s α κ c α i,κ ,
for all κ ∈ D \ I, and by

n i,Kσ = α∈P ζ α (p α Kσ , c α Kσ ) (φ K,Kσ s α K,Kσ + φ σ,Kσ s α σ,Kσ ) c α i,Kσ ,
for all interfaces Kσ ∈ I.

Discrete fluxes

The component discrete fluxes

q i,κ,κ = α∈P V α i,κ,κ + F α i,κ,κ ,
are the sum over each phase of the Darcy and Fickian diffusion fluxes which are defined as in ( 10)-( 12) by two-point flux approximations of the normal derivatives combined with a phase potential upwinding of the mobilities for the Darcy fluxes and an harmonic averaging of the saturations for the Fickian diffusion fluxes.

For convenience, for a given physical law ξ α function typically of p α and c α , we use the short notation ξ α κ = ξ(p α κ , c α κ ). The arithmetic mean will be denoted for brevity by

{{ξ α }} κ,κ = 1 2 (ξ α κ + ξ α κ ).
For ε > 0, we will also use a ε-regularized harmonic mean of two reals a ≥ 0 and b ≥ 0 given by the following expression

H ε (a, b) = 2ab a + b + ε . ( 20
)
Note that H 0 is simply the usual harmonic mean for a + b = 0.

Matrix-matrix fluxes -Let K, L ∈ M be two neighbouring cells sharing the face σ with σ ∈ F Γ . Let us define the transmissivities

T KL = |σ|H 0 1 x K -xσ , 1 x L -xσ and T D KL = |σ|H 0 ( Λ K x K -xσ , Λ L x L -xσ ),
then, we define the Darcy fluxes by

V α i,K,σ = ζ α K k α r,m (s α K ) µ α K c α i,K (V α K,σ ) + -ζ α L k α r,m (s α L ) µ α L c α i,L (V α K,σ ) - where V α K,σ = T D KL (p α K -p α L + {{ρ α }} K,L g(z K -z L )
) and the Fickian diffusion fluxes by

F α i,K,σ = D α T 2 m {{φ}} K,L {{ζ α }} K,L H ε (s α K , s α L )T KL (c α i,K -c α i,L ).
Fracture-fracture fluxes -Let σ, σ ∈ F Γ be two neighbouring fracture faces sharing the edge e ∈ E Γ \E . Let us define the transmissivities

T σσ = |e|H 0 1 xσ-xe , 1 x σ -xe and T D σσ = |e|H 0 ( Λσ xσ-xe , Λ σ x σ -xe ),
then, we define the Darcy fluxes by

V α i,σ,e = ζ α σ k α r,f (s α σ ) µ α σ c α i,σ (V α σ,e ) + -ζ α σ k α r,f (s α σ ) µ α σ c α i,σ (V α σ,e ) - where V α σ,e = T D σσ (p α σ -p α σ + {{ρ α }} σ,σ g(z σ -z σ ))
, and the Fickian diffusion fluxes by

F α i,σ,e = D α T 2 f {{φ}} σ,σ {{ζ α }} σ,σ H ε (s α σ , s α σ )T σσ (c α i,σ -c α i,σ ).
V α i,σ,e = ζ α σ k α r,f (s α σ ) µ α σ c α i,σ (V α σ,e ) + -ζ α e k α r,f (s α e ) µ α e c α i,e (V α σ,e ) - where V α σ,e = T D σe (p α σ -p α e + {{ρ α }} σ,e g(z σ -z e ))
, and the Fickian diffusion fluxes by

F α i,σ,e = D α T 2 f {{φ}} σ,e {{ζ α }} σ,e H ε (s α σ , s α e )T σe (c α i,σ -c α i,e ).
Matrix fracture interface fluxes -Let Kσ ∈ I with K ∈ M σ and σ ∈ F Γ . Let us define the half transmissivities

T K,Kσ = |σ| x K -xσ , T D K,Kσ = Λ K T K,Kσ , and
T σ,Kσ = 2 |σ| d f (xσ) , T D σ,Kσ = λ f,n,σ T σ,Kσ
, Then, we define the Darcy fluxes on the matrix and fracture sides by

V α i,K,Kσ = ζ α K k α r,m (s α K ) µ α K c α i,K (V α K,Kσ ) + -ζ α Kσ k α r,m (s α K,Kσ ) µ α Kσ c α i,Kσ (V α K,Kσ ) - and V α i,σ,Kσ = ζ α σ k α r,f (s α σ ) µ α σ c α i,σ (V α σ,Kσ ) + -ζ α Kσ k α r,f (s α σ,Kσ ) µ α Kσ c α i,Kσ (V α σ,Kσ ) - where V α K,Kσ = T D K,Kσ (p α K -p α Kσ + {{ρ α }} K,Kσ g(z K -z σ )), and V α σ,Kσ = T D σ,Kσ (p α σ -p α Kσ ) + {{ρ α }} σ,Kσ |σ|λ f,n,σ g • n Kσ .
The Fickian diffusion fluxes are defined on the matrix and fracture sides by

F α i,K,Kσ = D α T 2 m φ K {{ζ α }} K,Kσ H ε (s α K , s α K,Kσ )T K,Kσ (c α i,K -c α i,Kσ )
and

F α i,σ,Kσ = D α T 2 f φ σ {{ζ α }} σ,Kσ H ε (s α σ , s α σ,Kσ )T σ,Kσ (c α i,σ -c α i,Kσ ).

Numerical experiments

The first objective of this numerical section is to validate in Subsections 4.1 and 4.2 the hybrid-dimensional compositional two-phase Darcy flow model. This validation is performed by comparison with the reference equi-dimensional model obtained by a volumetric representation of the fractures. Two dimensional Cartesian DFM geometries and meshes are used for these validation test cases in order to be able to refine the mesh in the width of the fractures even for small apertures. Note that the TPFA discretization of the equi-dimensional model also incorporates face unknowns at mf interfaces to provide a physically consistent approximation of the mf interfacial fluxes.

Our hybrid-dimensional model is also compared with a classical reduced model obtained by a simple harmonic averaging of the half transmissivities at mf interfaces combined with an upwind approximation of the mobilities between the neighbouring cell and fracture face unknowns [START_REF] Karimi-Fard | An efficient discrete-fracture model applicable for general-purpose reservoir simulators[END_REF][START_REF] Gläser | A discrete fracture model for twophase flow in fractured porous media[END_REF][START_REF] Gläser | A hybrid-dimensional discrete fracture model for non-isothermal two-phase flow in fractured porous media[END_REF]. According to the terminology developped in [START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF] for immiscible two-phase DFM models, this harmonic averaging approximation will be termed in the following as mf linear f-upwind hybrid-dimensional model as opposed to the mf nonlinear hybrid-dimensional model considered in this work based on a nonlinear approximation at mf interfaces.

To be more specific, this model is derived directly at the discrete level using mf interface fluxes jumping over the interface. Let us define the following harmonic averages of the half transmissivities

T σ,K = H 0 ( φ σ T 2 f T σ,Kσ , φ K T 2 m T K,Kσ ), T D σ,K = H 0 (T D σ,Kσ , T D K,Kσ ),
then, these fluxes are defined by

V α i,σ,K = ζ α σ k α r,f (s α σ ) µ α σ c α i,σ (V α σ,K ) + -ζ α K k α r,m (s α K ) µ α K c α i,K (V α σ,K ) -,
for the Darcy fluxes, where

V α σ,K = T D σ,K (p α σ -p α K + {{ρ α }} σ,K g(z σ -z K )) ,
and by

F α i,σ,K = D α {{ζ α }} σ,K H ε (s α σ , s α K )T σ,K (c α i,σ -c α i,K
). for the Fickian fluxes. It is possible to derive the corresponding continuous mf interface fluxes by passing to the limit in the mesh size on the matrix side, leading to the mf interface Fickian fluxes

F α,± i,f m = ζ α (p α f , c α f ) + ζ α (γ ± p α m , γ ± c α m ) 2 s α f S α m (τ ± ) s α f + S α m (τ ± ) φ f D α T 2 f (c α i,f -γ ± c α i,m ) d f /2 ,
and to the mf interface Darcy fluxes

V α,± i,f m = c α i,f ζ α µ α (p α f , c α f ) k α r,f (s α f ) (V α,± f m ) + -γ ± c α i,m ζ α µ α (γ ± p α m , γ ± c α m ) k α r,m (S α m (τ ± )) (V α,± f m ) -, with V α,± f m = λ f,n (p α f -γ ± p α m ) d f /2 + ρ α (p α f , c α f ) + ρ α (γ ± p α m , γ ± c α m ) 2 g f • n ± .
By comparison with our reduced model based on the mf interface fluxes ( 10)-( 12), these limit fluxes exhibit a physical inconsistency in the sense that they use the matrix side interface saturation S α m (τ ± ) and relative permeability k α r,m rather than the fracture side saturation S α f (τ ± ) and relative permeability k α r,f . Let us however point out that practical meshes will never reach this limit due to the small fracture aperture and the usual high fracture permeability.

Two test cases are considered with liquid and gas phases defined as a mixture of air and water components. The first one presented in Subsection 4.1 simulates the gas injection from the bottom fracture in an initially liquid pure water saturated fractured porous medium. The fractures act as highly permeable drains. The gas rises mainly by gravity with strong capillary barrier effect at mf interfaces. Depending on the Henry coefficient, the dissolution of the air component in the liquid phase has a strong influence on the extension of the gas phase in the matrix domain while the Fickian diffusion is too low in the liquid phase to have a significant impact. The second test case of Subsection 4.2 investigates the effect of the Fickian diffusion in the gas phase at the mf interfaces. It considers the desaturation by suction of a fractured low permeable porous medium at the interface with the atmosphere.

A second objective of this numerical section is to consider in Subsection 4.3 the application of the hybrid-dimensional model to the simulation of the desaturation by suction of a low permeable and fractured storage rock. The set up of this test case is based on Andra nuclear waste storage prototype facility in a Callovo-Oxfordian argilite. The DFM is based on a simplified 2D version of the fracture network taking place at the interface between the exacavated tunnels and the storage rock.

In all the following test cases we consider a liquid wetting phase denoted by ( ) and a gas non-wetting phase denoted by (g). The liquid and gas phases are a mixture of two components i ∈ C = {air, h 2 o} including the water denoted by h 2 o which can vaporize in the gas phase, and the air component denoted by air which can dissolve in the liquid phase. The molar masses of the water and air components are fixed to m air = 29 • 10 -3 kg. mol -1 and m h2o = 18 • 10 -3 kg. mol -1 respectively. The gas molar density is defined by the perfect gas law ζ g = p g RT , with R = 8.314 J. K -1 . mol -1 and the liquid molar density is fixed to ζ = 1000 0.018 mol. m -3 . The phases viscosities are fixed to µ g = 18.51 • 10 -6 Pa. s -1 and µ = 10 -3 Pa. s -1 . The gas fugacities are given by Dalton's law for an ideal mixture of perfect gas f g i = c g i p g , i ∈ C. The fugacity of the air component in the liquid phase is given by Henry's law f air = c air H air (T ) with the temperature dependent Henry coefficient H air (T ) in Pa. For the water component in the liquid phase, the fugacity is taken from [START_REF] Schmidt | Properties of water and steam in S.I. units[END_REF]

f h2o = c h2o p sat (T ) exp - p sat (T ) -p ζ l RT ,
where p sat (T ) is the vapor pressure of the pure water given in Pa by the Rankine formula p sat (T ) = 1.013 • 10 5 exp 13.7 -5120 T .

We assume that the temperature is constant and fixed to T = 300 K, and the harmonic mean parameter ε in ( 20) is set to ε = 0.01. For all test cases, the parametrization of the capillary pressure graphs used at mf interfaces is based on the switch of variables between the fracture and matrix gas saturations (see for example Figure 2).

The simulation time is denoted by T f with an adaptive time stepping defined by ∆t 1 = ∆t init and for all n ≥ 1 by

∆t n+1 = max(∆t max , 1.2∆t n ),
in case of a successful time step, and ∆t n+1 = ∆t n 2 , in case of non convergence of the Newton algorithm in N L max = 50 iterations. The Newton solver is convergent if the relative residual is lower than 10 -5 or if the weighted maximum norm of the Newton increment is lower than 10 -4 . To obtain a more robust convergence of the nonlinear solver, a global damping of the Newton step is applied enforcing a prescribed maximum variation of the saturation. The Newton update is also locally modified at mf interfaces to avoid the jump in a single iteration across the point τ ± = 1 characterized by the discontinuity of dPc dτ ± , dS g ν dτ ± , ν ∈ {m, f }. These strategies are applied for all test cases.

At each Newton iteration, the linear system is solved using SuperLU 4.3 as a direct solver. The tests have been run sequentially on a laptop powered by an Intel Core i7 CPU clocked at 2.6 GHz and equipped with 32 Gb of RAM. We denote by N ∆t the number of successful time steps, by N Chop the number of time step chops, by N N ewton the average number of Newton iterations per successful time step. Finally, Cpu Time [s] stands for the CPU time in seconds.

Gas injection in porous fractures

We consider in this test case the DFM model exhibited in Figure 5 on the domain Ω = (0, 100 m) 2 . The Cartesian mesh used for the hybrid-dimensional models is uniform of size The matrix and fracture relative permeabilities and capillary pressures are given by k α r,f (s α ) = s α , k α r,m (s α ) = (s α ) 2 , and In order to investigate the effect of the dissolution on the gas saturation propagation, two values of the Henry coefficient are considered providing a high dissolution for H air = 10 8 Pa or a low dissolution for H air = 10 9 Pa. Figures 6 and7 obtained by the reference equi-dimensional model and the mf nonlinear and linear hybriddimensional models. It is clearly seen that the mf nonlinear model fully matches with the equi-dimensional model while differences appear with the mf linear model due to the rough approximation at mf interfaces. Basically, the mf linear f-upwind model underestimates the barrier effect leading to a higher leakage in the matrix as seen in Figure 6 and a slower gas propagation in the fractures combined with a higher air dissolution as exhibited in Figure 7.

P c,ν (s g ) = -b ν log 1 -s g (1 -e -aν bν ) , (21)
The numerical performances exhibited in Table 1 show that the mf nonlinear hybriddimensional model is as expected more efficient than the equi-dimensional model. The mf linear f-upwind model is penalized by many time step chops for the high dissolution case.

Desaturation by suction

We consider in this test case the DFM model exhibited in Figure 8 on the domain Ω = (0, 100 m) 2 . As for the previous test case, the Cartesian mesh is uniform of size n x × n y

The matrix and fracture domains are initially saturated by a pure water liquid phase. The top boundary conditions at y = 100 m are defined by pure water, the liquid saturation s m = 1 and the liquid pressure p m = 4 • 10 6 Pa.

The bottom boundary conditions at y = 0 are defined by p g m = p g f = 10 5 Pa and the relative humidity H r = 0.5 from which we deduce from the thermodynamical equilibrium the following boundary conditions both at the bottom matrix and fracture boundaries 9a and10, the mf linear f-upwind model does not capture this barrier removal effect due to a rough approximation of the matrix fracture interfacial Fickian diffusion fluxes. It is striking to see that the mf linear f-upwind model provides the same solution with and without Fickian diffusion. This is expected since, for the mf linear f-upwind model no diffusion can occur in the gas phase between the fracture and the cell above the fracture as long as the gas does not appear in the cell, which does not occur due to the barrier effect of gas filled fractures. On the other hand, the mf nonlinear model is able to model the gas diffusion flux between the fracture and the upper interface and to account for the appearance of the gas phase at the interface which then is transported to the upper cell by the gas Darcy flux (see paragraph 4.2.1 for more details). Table 2 shows the numerical performance of the three models with and without Fickian diffusion. It appears that the mf nonlinear hybrid-dimensional model is much more efficient than the equi-dimensional model thanks to a lower number of time step chops. It is also, as could be expected, somewhat more expensive in terms of Newton iterations than the mf linear f-upwind model due to its additional nonlinearities at mf interfaces. The data set is the same as in the previous test case except that we consider three different configurations w.r.t. Fickian diffusion.

       c g h2o = Hrpsat(T ) p g , c g air = 1 -c g h2o , c air = p g c g air Hair(T ) , c h2o = 1 -c air , p = p sat (T ) -ζ RT log( c h 2 o Hr ), s g ν = P -1 c,ν (p g -p ), ν ∈ {m, f }, (23)
• (nofick) : the Fickian diffusion is set to zero for both phases, i.e. D α = 0 for α ∈ P,

• (partialfick) : the Fickian diffusion is set to zero except for the gas phase and in the fractures only by setting T m → ∞ and D g = 2.5 • 10 -5 m 2 .s -1 ,

• (fick) : the Fickian diffusion is set to D g = 2.5 • 10 -5 m 2 .s -1 and D = 3 • 10 -9 m 2 .s -1 both in the matrix and in the fractures.

These three configurations are run using the mf nonlinear hybrid-dimensional model which provides basically the same solution as for the equi-dimensional model. Figures 12 and13 exhibit that the solutions obtained by the (fick) and (partialfick) cases fully match showing that the gas diffusion in the fracture is responsible for the removal of the barrier effect clearly observed for the (nofick) case. This can be understood from Figure 14 which exhibits, as a function of time, the parameter τ at the upper left horizontal interface, the Darcy flux for the air component in the gas phase between this interface and the upper cell, and the Fickian flux for the air component in the gas phase between the fracture and this interface. It shows that the air component is transported by gas diffusion from the fracture to the upper interface which triggers the appearance of the gas phase at the interface on the matrix side (parameter τ > 1). This gas phase is then transported by the gas Darcy flux to the upper cell, see Figure 14a for an illustration with the convention that a positive flux is directed towards the interface.

Application to the desaturation by suction of a Callovo-Oxfordian argilite fractured rock

The rectangular domain Ω = (0, 5 m) × (0, 10 m) and the fracture network of this application test case are illustrated in Figure 15 together with the triangular mesh generated with Gmsh [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with builtin pre-and post-processing facilities[END_REF]. The bottom boundary represents the interface between the matrix and the ventilation tunnel with a relative humidity fixed to H r = 0.6. Within around one meter from the gallery, the rock has been fractured by the excavation of the tunnel. The resulting Excavation Damaged Zone (EDZ) presents oblique "chevron" fractures and unloading joints 12 and the porosity set to φ f = 1. The matrix relative permeabilities and capillary pressure are given by the following Van Genuchten laws with the parameters n = 1.49, m = 1 -1 n , P r = 15 • 10 6 Pa and the residual liquid and gas saturations s r = 0.4, s g r = 0:

k r,m (s ) =        0 if s < s r , 1 if s > 1 -s g r , √ s 1 -(1 -(s ) 1 m ) m 2 if s r ≤ s ≤ 1 -s g r , V g air,K,Kσ ≤ 0 F g air,σ,Kσ ≥ 0 K σ (
g r,m (s g ) =        0 if s g < s g r , 1 if s g > 1 -s r , √ 1 -s 1 -(s ) 1 m 2m if s g r ≤ s g ≤ 1 -s r , P c,m (s ) = P r ((s ) -1 m -1) 1 n if 0 < s ≤ 1, with s = s -s r 1 -s r -s g r .
The fracture relative permeabilities are defined by k α r,f (s α ) = s α and the fracture capillary pressure is given by ( 21) with b f = 10 2 Pa, a f = 10b f .

The Fickian diffusion coefficients are set to either D g = D = 0 m 2 .s -1 or to D g = 2.5 • 10 -5 m 2 .s -1 and D = 3 • 10 -9 m 2 .s -1 in order to investigate the effect of the diffusion on the gas propagation. The porous medium tortuosities are fixed to T m = 2, T f = 1, and the Henry coefficient if given by [START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal darcy flows in fractured porous media[END_REF]. The matrix and fracture domains are initially saturated by a pure water liquid phase. The top boundary conditions are defined by pure water, the liquid saturation s m = 1 and the liquid pressure p m = 4 • 10 6 Pa. The bottom boundary conditions at y = 0 are defined both at the matrix and fracture boundaries by p g m = p g f = 10 5 Pa and the relative humidity H r is equal to 0.6 from which we deduce from the thermodynamical equilibrium the boundary conditions given by [START_REF] Brenner | Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF]. The lateral boundaries are assumed impervious using periodicity arguments. The simulation time is set to T f = 200 years and the time stepping is defined by the parameters ∆t 1 = 10 -5 s and ∆t max = 2 years. This test case is run with or without Fickian diffusion using the mf nonlinear hybrid-dimensional model.

Sensitivity of the solution to the normal fracture permeability

The normal fracture permeability is not properly defined for open fractures which raises the issue of how to choose this model parameter. Fortunately we check in this paragraph that the solution is only weakly sensitive to the value of the normal fracture permeability as long as it is much larger than the matrix permeability. Let us set

λ f,n = rΛ f . (24)
In order to investigate the sensitivity of the solution to the choice of the parameter r, we consider a single oblique fracture with aperture size d f = 10 mm. Figure 17 compares the gas volume, the water vapor and dissolved air number of moles both in the matrix and in the fracture as a function of time. This comparison is performed for three different values of the parameter r = 1, 10 -3 , 10 -6 , where r = 1 corresponds to an isotropic fracture permeability. In addition, Figure 16 plots the gas saturation for each value of r at final time T f = 200 years with and without Fickian diffusion. Small differences are observed on the fracture gas volume during the fracture infill at short times, then following the fracture infill, no significant differences are observed in the fracture nor in the matrix. The same remark applies to the number of moles of water vapor and of dissolved air. In overall, it appears that the impact of the normal fracture permeability is quite small for a large range of the parameter r, roughly speaking as long as the normal fracture permeability remains much larger than the matrix permeability. Table 3 shows the numerical performances for all these test cases. We observe as expected that smaller values of the parameter r are advantageous since it reduces the nonlinearity and hence the number of Newton iterations or time step chops especially for the no Fickian diffusion case, without affecting significantly the solution.

Full fracture network

We now consider the fracture network exhibited in the right Figure 15 with d f = 1 mm for the oblique fractures and d f = 0.01 mm for the horizontal fractures. The parameter r of the fracture normal permeability in ( 24) is fixed to r = 10 -6 . This choice is in the range of low sensitivity of the solutions to the normal fracture permeability and provide a good nonlinear convergence behavior. Figure 18 exhibits the gas saturation, the gas pressure, and the water vapor molar fraction at final time with and without the Fickian diffusion. Figure 19 plots the solutions over the vertical segment at x = 2.7 m with y ∈ (0, 5 m) crossing the EDZ (see the vertical green line in Figure 15b). The barrier effect induced by the gas filled fractures is highly visible on the saturation plots. As already noticed in the previous test case, this barrier effect is considerably reduced by the Fickian diffusion of the gas phase in the fractures. We also note that the gas Fickian diffusion homogenize the water vapor molar fraction and consequently the gas pressure in the gas region.

Without Fick

Table 3 shows in the Full network lines, the numerical performances obtained with or without diffusion where we can observe a higher number of time step chops for the diffusion case.

Conclusion

This article presents a new hybrid-dimensional compositional two-phase Darcy flow model accounting for phase transitions and Fickian diffusion. The nonlinear transmission conditions are designed to be consistent with the main physical processes at mf interfaces. They are based on the flux continuity of each component using Two-Point Flux Approximations (TPFA) in the width of the fracture. The saturation jumps are captured using a parametrization of the capillary pressure graphs. The thermodynamical equilibrium is formulated using complementary constraints and taking into account the saturation jumps at mf interfaces. The hybrid-dimensional model is discretized using a TPFA combined with a redistribution of the porous volume both at mf interfaces and at edges shared by more than three fracture faces. The hybrid-dimensional model is shown to provide basically the same accuracy as the reference equi-dimensional model and to be more accurate and physically consistent than the usual approach based on harmonic averaging of the transmissivities at mf interfaces combined with a two-point upwinding of the mobilities jumping over the mf interfaces. This hybrid-dimensional model is used to investigate the desaturation of a fractured Callovo-Oxfordian argilite at the interface with a ventilation tunnel. The Fickian diffusion in the gas phase and in the fracture width is shown to remove the barrier effect induced by the gas filled fractures. In perspective, the discretization will be extended to account for general meshes and anisotropy of the permeability tensor using a face based approach in the spirit of [START_REF] Aghili | Hybrid finite volume discretization of two-phase discrete fracture matrix models with nonlinear interface solver[END_REF].

Figure 1 :

 1 Figure 1: Example of a 2D domain Ω with three intersecting fractures Γ i , i = 1, 2, 3.

Figure 2 :

 2 Figure 2: Example of parametrization of the capillary pressure graphs based on a switch of variables between the fracture and matrix non-wetting phase saturations. Plots of S nw ν (τ ), ν ∈ {m, f } and of Pc(τ ) am for the capillary pressures P c,ν (s g ) = -b ν log 1 -s g (1 -e -aν bν ) with b f = 10 5 Pa, b m = 10 6 Pa, a f = 2b f , a m = 5b m , τ ∈ (τ 0 , τ 1 ) with τ 0 = 0, τ 1 =

Figure 4 :

 4 Figure 4: Discrete unknowns of the compositional two-phase flow model on an orthogonal triangular mesh discretized with the two-point flux approximation. The fracture faces are in bold. The mesh orthogonality assumption requires acute angles and to choose the cell centre x K at the intersection of the orthogonal bisectors, and the face centre x σ at the face middle point.

Fracture-edge

  fluxes -Let σ ∈ F Γ and e ∈ E σ ∩E . Let us define the half transmissivities T σe = |e| xσ-xe and T D σe = Λ σ T σe , then, we define the Darcy fluxes by

Figure 5 :

 5 Figure 5: Gas injection test case -Discrete Fracture Model with the fracture network in bold and its uniform Cartesian mesh of size 16 × 16.

with b f = 5 • 10 3

 510 Pa, b m = 10 6 Pa, a f = 2b f , a m = 5b m . The matrix is homogeneous and charaterized by the isotropic permeability Λ m = 0.01 Darcy, and the porosity φ m = 0.2. The fracture zone properties are set to d f = 1 cm, Λ f = λ f,n = 10 4 Λ m , φ f = 0.4. The Fickian diffusion coefficients are set to D g = 2.5 • 10 -5 m 2 .s -1 and D = 3 • 10 -9 m 2 .s -1 , and the porous medium tortuosities to T m = T f = 1. The matrix and fracture domains are initially saturated by a pure water liquid phase. The top boundary conditions at y = 100 m are defined by pure water, the liquid saturation s m = 1 and the liquid pressure p m = 4 • 10 6 Pa. The bottom boundary conditions at y = 0 are impervious in the matrix and defined in the bottom fracture by the water pressure p f = 4 • 10 6 + 100g ρ Pa, with ρ = 1000 kg. m -3 , the gas saturation s g f = 1 and pure air composition. The lateral boundaries are assumed impervious.

  The simulation time is set to T f = 100 days. The time stepping is defined by the parameters ∆t 1 = 1 s, ∆t max = 0.05 day for t ≤ 10 days and ∆t max = 1 day for t > 10 days. H air = 10 8 (a) Equi-dim (b) mf nonlinear (c) mf linear f -upwind H air = 10 9 (d) Equi-dim (e) mf nonlinear (f) mf linear f -upwind

Figure 6 :

 6 Figure 6: Gas injection test case -Comparison of the gas saturation at final time using the three different models: equi-dimensional model (left), hybrid-dimensional mf nonlinear model (center) and hybrid-dimensional mf linear f -upwind model (right) and the two values of the Henry coefficient H air .

Figure 7 :

 7 Figure 7: Gas injection test case -For the three models and the two Henry coefficients, gas volume (left), water vapor number of mole (center), and dissoluted air number of mole (right) as a function of time (in days) both in the matrix (top) and in the fractures (bottom).

Figure 8 :

 8 Figure 8: Discrete Fracture Model with the fracture network in bold and its uniform Cartesian mesh of size 16 × 16.

resulting in roughly s g m = 0.86 and s g f = 1 .

 1 The lateral boundaries are assumed impervious. The simulation time is set to T f = 200 years and the time stepping is defined by the parameters ∆t 1 = 10 -3 s and ∆t max = 2 years.

  (a) mf nonlinear hybrid-dim vs equi-dim (b) mf nonlinear vs. mf linear f-upwind

Figure 9 :Figure 10 :

 910 Figure 9: Suction test case -For the three models and with or without Fickian diffusion, gas volume (left), water vapor number of mole (center), and dissoluted air number of mole (right) as a function of time (in days) both in the matrix or in the fractures.

4. 2 . 1

 21 Removal by Fickian diffusion of the barrier effect induced by gas filled fracturesIn order to better understand the role of the Fickian diffusion on the removal of the barrier effect induced on the liquid phase by gas filled fractures, let us consider the simplified DFM model exhibited in Figure11on the domain Ω = (0, 20 m) × (0, 100 m) meshed by a uniform 2 × 128 Cartesian mesh.

Figure 11 :

 11 Figure 11: Discrete Fracture Model with the fracture network in bold and its uniform Cartesian mesh of size 2 × 128.

Figure 12 :

 12 Figure 12: Two fractures suction test case -Gas saturation fields at final time for the (nofick) case (left), the (partialfick) case (center), and the (fick) case (right).

Figure 13 :

 13 Figure 13: Two fractures suction test case -Matrix and fracture gas volume as a function of time for the (nofick), (partialfick) and (fick) cases.

  a) Fluxes at the upper interface Kσ between a cell K and a fracture face σ. For readability, the fracture is represented with the same dimension as K and with an enlarged aperture. (b) Variable tau at upper interface (c) Darcy flux V g air,K,Kσ (d) Fickian flux F g air,σ,Kσ

Figure 14 :

 14 Figure 14: Two fractures suction case -For the (nofick), (partialfick) and (fick) cases: the parameter τ at the upper left horizontal interface Kσ ( ), the Darcy flux V g air,K,Kσ for the air component in the gas phase between the interface Kσ and the upper cell K (represented by ♦), and the Fickian flux F g air,σ,Kσ for the air component in the gas phase between the fracture σ ( ) and the interface Kσ.

  k

C

  = {air,water} P = { , g} (a) Domain Ω with the EDZ (b) Computational mesh

Figure 15 :

 15 Figure 15: Discrete Fracture Matrix model with the fracture network in the damaged zone (EDZ) and its triangular mesh. The boundary condition at the interface between the porous medium and the ventilation tunnel is defined by the atmospheric pressure p g = 10 5 Pa and the relative humidity H r = 0.6.

(a) r = 1

 1 (b) r = 10 -3 (c) r = 10 -6 With Fick (d) r = 1 (e) r = 10 -3 (f) r = 10 -6

Figure 16 :

 16 Figure 16: Comparison on the one oblique fracture test case of the gas saturation fields at final time for different values of the normal fracture permeability obtained with r = 1, 10 -3 , 10 -6 and with (bottom) and without (top) Fickian diffusion.

Figure 17 :

 17 Figure 17: Comparison on the one oblique fracture test case of the gas volume, and water vapor number of moles, and dissolved air number of moles for r = 1, 10 -3 , 10 -6 and the cases with or without Fickian diffusion.

( a )

 a Gas saturation where s g > 10 -6 (b) Gas pressure where 0.9 • 10 5 ≤ p g ≤ 1.1 • 10 5 Pa (c) Water vapor molar fraction c g h 2 o , zoom in the EDZ

Figure 18 :

 18 Figure 18: Andra test case -Zoom on regions of the gas saturation, pressure and water vapor molar fraction at final time with Fickian diffusion (left) and without Fickian diffusion (right).

Figure 19 :

 19 Figure 19: Andra test case -Cuts of the saturation, pressure and water vapor at final time (see Fig. 18) along a vertical line crossing the EDZ from (2.7, 0) to (2.7, 5), see Fig.15b.

Table 1 :

 1 compare the solutions Test case H air Nb Cells Cpu Time [s] N ∆t N Chop N N ewton Gas injection test case -Numerical performances, Nb Cells is the number of mesh cells, CPU time is the total execution time in seconds, N ∆t is the number of successful time steps, N N ewton is the average number of Newton iterations per successful time step, and N Chop is the number of time step failures.

	mf nonlinear	10 8	4096	99.25	342	0	5.47
	mf nonlinear	10 9	4096	127.5	344	1	5.01
	mf linear f -upwind 10 8	4356	283.0	574	75	9.46
	mf linear f -upwind 10 9	4356	65.1	342	0	3.45
	Equi-dim	10 8	5476	276.64	380	16	10.21
	Equi-dim	10 9	5476	245.94	368	12	9.23

Table 2 :

 2 Test caseFick diff. Nb Cells Cpu Time [s] N ∆t N Chop N N ewton mf linear f -upwind Suction test case -Numerical performances, Nb Cells is the number of mesh cells, CPU time is the total execution time in seconds, N ∆t is the number of successful time steps, N N ewton is the average number of Newton iterations per successful time step and N Chop is the number of time step failures. boundary. The Fickian diffusion considerably reduces this barrier effect. This physical process is investigated in more details in the following paragraph 4.2.1. As shown in Figures

		with	4356	72.56	329	11	4.33
	mf linear f -upwind without	4356	119.84	373	26	6.09
	mf nonlinear	with	4096	179.45	414	21	8.02
	mf nonlinear	without	4096	161.97	357	9	8.4
	Equi-dim	with	5476	351.07	443	34	10.73
	Equi-dim	without	5476	623.5	530	62	15.83

Table 3 :

 3 Andra test case -Numerical performances, Nb Cells is the number of mesh cells, CPU time is the total execution time in seconds, N ∆t is the number of successful time steps, N N ewton is the average number of Newton iterations per successful time step and N Chop is the number time step failures.

	Full network	with	10 -6	11549	403.08	384	34	7.18
	Full network	without 10 -6	11549	259.6	280	6	5.62
	One fracture	with	1	11457	138.83	257	0	3.7
	One fracture without	1	11457	415.9	325	18	9.02
	One fracture	with	10 -3	11457	120.11	259	1	3.15
	One fracture without 10 -3	11457	349.61	318	16	8.08
	One fracture	with	10 -6	11457	109.8	259	1	2.82
	One fracture without 10 -6	11457	103.32	257	0	2.49

Test case

Fick diff.

r Nb Cells Cpu Time [s] N ∆t N Chop N N ewton

Acknowledgement:

The authors would like to thank Andra, the french national radioactive waste management agency, for supporting this work.


Pierre Simonin

Carine Lombard

Arnaud Huguet

Adrienne Kish
email: adrienne.kish@mnhn.fr

Improved Isolation of SlaA and SlaB S-layer Proteins in Sulfolobus acidocaldarius

Keywords: S-layer, cell wall, Sulfolobus acidocaldarius, Archaea, protein purification, SlaA, SlaB

SDS = sodium dodecyl sulfate SDS-PAGE = sodium dodecyl sulfate polyacrylamide gel electrophoresis HEPES = 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid EDTA = Ethylenediaminetetraacetic acid MeOH = Methanol DCM = Dichloromethane HPLC-APCI-MS = high performance liquid chromatography

The hyperthermoacidophilic Sulfolobus acidocaldarius has served as a key model for studying the S-layers of archaea since its isolation in 1972 [START_REF] Brock | Sulfolobus: a new genus of sulfuroxidizing bacteria living at low pH and high temperature[END_REF]. Like other Sulfolobus species, the S-layers of S. acidocaldarius are composed of monomers of the SlaA protein self-assembled into a hexagonal array with threefold (P3) symmetry [START_REF] Lembcke | Image analysis and processing of an imperfect two-dimensional crystal: the surface layer of the archaebacterium Sulfolobus acidocaldarius re-investigated[END_REF], and anchored into a tetraether monolayer cell membrane by the SlaB protein proteins [START_REF] Veith | Acidianus, Sulfolobus and Metallosphaera surface layers: structure, composition and gene expression[END_REF][START_REF] Gambelli | Architecture and modular assembly of Sulfolobus S-layers revealed by electron cryotomography[END_REF]. Both SlaA and SlaB are known to be Nglycosylated [START_REF] Peyfoon | The S-layer glycoprotein of the crenarchaeote Sulfolobus acidocaldarius is glycosylated at multiple sites with chitobiose-linked Nglycans[END_REF][START_REF] Guan | N-glycosylation in the thermoacidophilic archaeon Sulfolobus acidocaldarius involves a short dolichol pyrophosphate carrier[END_REF]. While Sulfolobus species share the same the basic S-layer structure, the primary sequence of SlaA and SlaB share only 24-25 % identify with other species such as S. islandicus and S. solfataricus {Gambelli:2019hd}. The S-layer proteins of S. acidocaldarius also show far greater stability against detergents than other Sulfolobus species. While Sulfolobus solfataricus and Sulfolobus brierleyi are only mildly resistant to hot detergent (SDS) treatments, S. acidocaldarius was shown to be highly resistant [START_REF] König | Studies on archaebacterial S-layers[END_REF]. Treatment with room temperature 10 % SDS resulted in Slayer sheets, but the S. acidocaldarius S-layer proteins remained insoluble up to treatment with 10 % SDS at 90 o C [START_REF] Weiss | Subunit cell wall of Sulfolobus acidocaldarius[END_REF]. Many aspects of the S-layer structure and function can be studied using whole cells, however isolation and purification of the S-layer proteins is needed to separate the roles of SlaA and SlaB from those of the other cell envelope components. Despite the simplistic structure of the Sulfolobus cell envelope, isolation of the S-layer proteins is complicated by the fact that after an initial cell disruption and removal of the cellular contents, S-layer-containing cell-free sacculi retain the same morphology as whole cells. These sacculi or "ghosts" can contain the complete cell envelope (cell membrane, SlaA, SlaB), or only the cell wall (S-layer only, no cell membrane). Early extractions of the S. acidocaldarius cell envelope lysed cells using either a French press, or detergents (0.15 % SDS with 10 mM EDTA) to aid in destabilizing the cell envelope [START_REF] Weiss | Subunit cell wall of Sulfolobus acidocaldarius[END_REF][START_REF] Michel | The 2-D Crystalline Cell Wall of Sulfolobus acidocaldarius: Structure, Solubilization, and Reassembly[END_REF]. [START_REF] Grogan | Phenotypic characterization of the archaebacterial genus Sulfolobus: comparison of five wild-type strains[END_REF] developed a generalized Sulfolobus S-layer purification protocol using 0.4 % N-lauryl sarcosine in buffer (pH 5) to lyse cells and DNase treatment to reduce the lysate viscosity, followed by repeated cycles of treatment with the same solution heated to 56 o C and centrifugation to remove cellular debris, eventually producing a whitish upper phase containing the S-layer proteins overlaying the brownish cell pellet. The higher resistance of S. acidocaldarius S-layer proteins to detergent treatments compared to other Sulfolobus species were exploited to permitted better separation of the S-layer proteins from the underlying cell membrane. Lipid removal was found to be more effective using SDS rather than Triton X-100. Using large batches (over 20 L) of cell cultures for extractions, it was found that following an initial cell lysis, repeated hot (60 o C) SDS treatments could be employed to separate S. acidocaldarius S-layers from other cellular proteins [START_REF] Michel | The 2-D Crystalline Cell Wall of Sulfolobus acidocaldarius: Structure, Solubilization, and Reassembly[END_REF]. SDS was subsequently removed by 6-8 washes in deionized H2O (dIH2O). This protocol became the standard for S. acidocaldarius S-layer extractions. However, the composition of the protein extracts (SlaA, SlaB) varied between labs. Using the same protocol and the same strain of S. acidocaldarius, one group reported obtaining S-layer sacculi composed of both SlaA and SlaB proteins [START_REF] Michel | The 2-D Crystalline Cell Wall of Sulfolobus acidocaldarius: Structure, Solubilization, and Reassembly[END_REF], while others reported obtaining pure SlaA "ghosts" retaining cell-type morphology [START_REF] Selenska-Pobel | Magnetic Au nanoparticles on archaeal S-Layer ghosts as templates[END_REF][START_REF] Reitz | Spectroscopic study on uranyl carboxylate complexes formed at the surface layer of Sulfolobus acidocaldarius[END_REF]. The quantity of cell culture used for extraction also affects the quality of the extracts. Following the same protocol, but using only 1 L of exponential-phase culture instead of large batch cultures (20-70 L), we obtained S-layer ghosts with both SlaA and B proteins rather than pure SlaA [START_REF] Kish | Preservation of Archaeal Surface Layer Structure During Mineralization[END_REF]).

P r e -P r i n t p r i o r t o p e e r r e v i e w

s e e f i n a l v e r s i o n d o i : 1 0 . 1 0 0 7 / s 0 0 7 9 2 -0 2 0 -0 1 1 7 9 -9

The inability to properly down-scaling the S-layer extraction protocol to benchtop scale (1-2 L) cultures, as well as the difficulties encountered isolating SlaB, likely due to the strong association of SlaB with so-far unidentified molecules in the cell envelope, inhibit routine studies of S-layer protein structure and function.

Here we present a method for a benchtop-scale purification of SlaA and SlaB proteins, as well as characterization of the molecules strongly associated with the SlaB anchor protein.

Methods and Results:

The purification of SlaA and SlaB proteins followed three main steps: preparation of cell envelope extracts, separation of SlaA and SlaB, and finally removal of membrane lipids from SlaB.

Culture growth and preparation of S-layer-enriched cell envelope extracts Sulfolobus acidocaldarius DSM 639 cells were cultivated in Brock's medium [START_REF] Brock | Sulfolobus: a new genus of sulfuroxidizing bacteria living at low pH and high temperature[END_REF]) pH 3.5 supplemented with 0.1% yeast extract and 0.2% D-saccharose at 80 o C with agitation (170 rpm) in 1 L Erlenmeyer flasks with specialty-made extended glass aeration necks in a shaking oil bath until mid-exponential growth phase (OD600nm=0.5). The cells were harvested by centrifugation (8000 x g, 10 min), and the cell pellets resulting from roughly 1 L of culture were flash frozen in liquid nitrogen and stored at -80 o C until use for S-layer protein extraction.

To purify the S-layer proteins, we first isolated the S-layer containing cell envelope using a modified version of the protocol of Reitz et al. [START_REF] Reitz | Spectroscopic study on uranyl carboxylate complexes formed at the surface layer of Sulfolobus acidocaldarius[END_REF], adapted to benchtop-scale volumes (1-2 L cell culture at OD600nm=0.5) rather than bioreactor-scale volumes (70 L) used in previous studies. The cell pellet was gently thawed on ice and resuspended in HEPES buffer (10 mM, pH 7) at 4 o C containing 2 mM EDTA to inhibit metalloproteases as well as chelate any cations attached to anionic leaflet of the tetraether membrane [START_REF] De Castro | Biological Membranes in Extreme Conditions: Anionic Tetraether Lipid Membranes and Their Interactions with Sodium and Potassium[END_REF], and 1 mM Pefabloc SC to inhibit serine proteases [START_REF] Fusi | Intracellular proteases from the extremely thermophilic archaebacterium Sulfolobus solfataricus[END_REF]. Preliminary disruption of cell membranes and solubilization of proteins was achieved using a denaturing, anionic detergent (0.15 % SDS). DNA was digested by addition of 100 µg/µl DNase I and 4 mM MgCl2 for 1 h at 37 o C to reduce the viscosity of the cell lysate. Further solubilization of the cell lysate and cell membrane was then performed by overnight incubation in 2 % SDS at room temperature with gentle agitation (50 rpm). Insoluble S-layer components were separated from the rest of the solubilized cell lysate by centrifugation (40 000 x g, 45 min). The white upper, diffuse layer containing the cell envelope was removed from the more compact, darker beige phase underneath, resuspended in the same buffer, and centrifuged as before. The SDS was then removed from the resulting cell layer extracts by eight washes in dIH2O (centrifugation at 40 000 x g, 45 min between each wash step), and the protein components of extracts analyzed by SDS-PAGE. This enriched S-layer protein sample obtained by hot SDS extraction was then analyzed by SDS-PAGE. Solubilization of the proteins prior to electrophoretic migration was accomplished by incubation in 20 mM sodium carbonate buffer (pH 10) at 65 o C for 30 min.

Laemmli buffer was then added and the proteins denatured at 95 o C for 5 min prior to separation by SDS-PAGE (10 % gel), and visualized by Coomassie staining. The only two bands visible (apparent masses near 120 kDa and 55 kDa) corresponded to the previously established apparent molecular masses for polypeptides determined by mass spectrometry to be SlaA and SlaB, respectively [START_REF] Kish | Preservation of Archaeal Surface Layer Structure During Mineralization[END_REF]) (see Figure 1). More SlaA was isolated than SlaB, which is coherent with the model for Sulfolobus S-layer structure, with six SlaA

P r e -P r i n t p r i o r t o p e e r r e v i e w

s e e f i n a l v e r s i o n d o i : 1 0 . 1 0 0 7 / s 0 0 7 9 2 -0 2 0 -0 1 1 7 9 -9

proteins composing the hexameric S-layer array anchored by three SlaB proteins [START_REF] Veith | Acidianus, Sulfolobus and Metallosphaera surface layers: structure, composition and gene expression[END_REF].

After removal of the supernatant, the resulting pellet contains two differently colored phases: a lower, darker beige phase containing cell lysate and a whitish upper phase containing the cell envelope. Previous protocols manually separated the whitish phase from the rest of the pellet [START_REF] Reitz | Spectroscopic study on uranyl carboxylate complexes formed at the surface layer of Sulfolobus acidocaldarius[END_REF]. However, in practice obtaining a cell pellet sufficiently large to easily manually separate the whitish phase from the darker phase is only possible when starting from batch-scale culture volumes, such as the 70 L bioreactor used by Reitz et al. [START_REF] Reitz | Spectroscopic study on uranyl carboxylate complexes formed at the surface layer of Sulfolobus acidocaldarius[END_REF], which are not readily available. To separate the SlaA and SlaB proteins, we instead used size exclusion chromatography.

Separation of SlaA and SlaB by size exclusion chromatography

The S-layer suspension in 20 mM Na-carbonate buffer (pH 10) was injected into a Superdex 200 (GE Healthcare) size exclusion chromatography column. Figure 2 shows the presence of two peaks: the first is an exclusion peak (panel A; peak 1), whereas the second one corresponds to purified SlaA monomers (panel A; peak 2). The purified SlaA monomers were collected. In order to ensure that no SlaA monomers were mixed with SlaB proteins in the excluded fraction, this fraction was reinjected and separated by gel filtration (Figure 2B; peak 1'). The resulting fractions from both injections were stored at 4 o C. Separation of the isolated proteins by SDS-PAGE (10 % gel; see Figure 3) revealed the expected apparent molecular masses for purified SlaB and SlaA proteins in the first and second fractions, respectively.

While this shows that the SlaA and SlaB proteins were successfully separated from all other proteins, the chromatogram for the elution of the SlaB protein corresponds to a larger mass than would be expected for SlaB monomers. We hypothesized that SlaB monomers were aggregating, either due to hydrophobic interactions between SlaB proteins, or due to the presence of membrane tetraether lipids still associated with the SlaB protein. To test the hypothesis that lipids were remaining in association with SlaB even after treatment with 2% SDS, we subjected the enriched SlaB fraction to lipid extraction. Extraction and identification of membrane lipids associated with SlaB The enriched SlaB suspension collected from the size exclusion chromatography column (Figure 2B, peak 1') was dried using a rotary evaporator. Lipid extraction from the dried sample was performed using a modified Bligh and Dyer method based on [START_REF] Huguet | Effects of a short-term experimental microclimate warming on the abundance and distribution of branched GDGTs in a French peatland[END_REF]. Briefly, samples were resuspended in an extraction mix of methanol (MeOH):Dichloromethane (DCM):Buffer (either phosphate buffer (50 mM, pH 7.4) or Nacarbonate buffer (40mM, pH 10)) (2:1:0.8 v/v/v), and placed in an ultrasonic bath (Branson 3510) for 10 min. The MeOH:DCM:Buffer (Phosphate or Na-carbonate) volume ratio was then adjusted to 1:1:0.9, followed by centrifugation (500 x g, 10 min), and collection of the extracts in a separation funnel. The DCM layer containing the lipids was separated from the MeOH: Buffer (Phosphate or Na-carbonate) layer and then extracted three more times with DCM, with the resulting DCM fractions combined and dried over Na2SO4. For extraction of lipids from the enriched SlaB suspension, Na-carbonate (40 mM, pH 10) was used.

The lipid extract was then re-dissolved in 200 µL heptane and centrifuged using an Eppendorf Mini Spin centrifuge (1 min, 7000 rpm). The supernatant was collected and analyzed using high performance liquid chromatography-atmospheric pressure chemical ionization mass spectrometry (HPLC-APCI-MS). Tetraether analyses were performed as

P r e -P r i n t p r i o r t o p e e r r e v i e w

s e e f i n a l v e r s i o n d o i : 1 0 . 1 0 0 7 / s 0 0 7 9 2 -0 2 0 -0 1 1 7 9 -9 recently described by Huguet et al. (2019). Briefly, tetraether separation was achieved with two silica columns in tandem (150 mm × 2.1 mm, 1.9 µm, Thermo Finnigan; USA) thermostated at 40 o C. Injection volume was 30 µl. GDGTs were first eluted isocratically with 82 % A/18 % B for 25 min (A = hexane, B = hexane/isopropanol 9/1, v/v). The following linear gradient was subsequently used: 82 % A/18 % B to 65 % A/35 % B in 25 min, followed by 65 % A/35 % B to 100 % B, maintained for 10 min and then back to 82 % A/18 % B in 10 min, maintained for 30 min. The flow rate was set at 0.2 ml/min. HPLC-APCI-MS analysis enabled detection of two tetraether lipids in the sample of purified SlaB proteins, identified on the chromatogram as two peaks with mass to charge ratios (m/z) of 1296 at 10.7 min, and 1294 at 11.4 min (Figure 4A). These masses correspond to the cyclopentane ring-bearing glycerol dialkyl glycerol tetraethers (GDGTs) GDGT-3 and GDGT-4, respectively [START_REF] Hopmans | Analysis of intact tetraether lipids in archaeal cell material and sediments by high performance liquid chromatography/atmospheric pressure chemical ionization mass spectrometry[END_REF][START_REF] Zeng | GDGT cyclization proteins identify the dominant archaeal sources of tetraether lipids in the ocean[END_REF].

To ensure that no other lipids remained attached to the SlaB proteins after lipid extraction, the SlaB proteins obtained after solvent-based lipid extraction were first applied to a Superdex 200 (GE Healthcare) size exclusion chromatography column. The resulting protein fraction was then subjected to a second round of solvent-based lipid extraction using the same protocol as previously, using phosphate buffer (50 mM, pH 7.4) instead of Nacarbonate buffer, after which no lipids were detected (Figure 4B) with a minimum detection limit of 1 ppb.

Having identified that GDGTs were in fact associated with SlaB after hot SDS treatment and size exclusion chromatography, we sought to remove these lipids from the original SlaB protein aggregate obtained after hot SDS treatment. This enriched SlaB protein solution (Figure 2B, peak 1') was first solubilized in an equal volume of 20 mM Na-carbonate buffer (pH 10) and then treated it with 1 % CHAPS for 30 min at 37 o C with agitation. To check whether this treatment was sufficient for removal of the GDGTs associated with SlaB, the CHAPS-treated solution was first applied to a Superdex 200 size exclusion chromatography column to isolate the SlaB proteins, which were excluded from the column, followed by lipid extraction following the same modified Bligh and Dyer method described above using phosphate buffer (50 mM, pH 7.4). HPLC-APCI-MS analyses failed to detected any tetraether lipids in association with SlaB. The exclusion of SlaB proteins from the Superdex 200 resin indicates that despite the removal of residual GDGTs by CHAPS following hot SDS isolation of the S-layer proteins, the SlaB proteins remained associated in an aggregate mass.

Taken together, these results suggest that the larger mass observed during elution for purified SlaB was due to both strong association with GDGTs and likely auto-aggregation of hydrophobic, N-terminal transmembrane domain of the SlaB protein [START_REF] Veith | Acidianus, Sulfolobus and Metallosphaera surface layers: structure, composition and gene expression[END_REF].

Discussion :

Here we present a method for complete isolation of SlaA from SlaB, starting from benchscale culture volumes. This method improves upon existing protocols by allowing for reduced quantities of cell cultures to be used. In our study, the SlaA was found to have an apparent molecular mass of 120 kDa, despite the calculated molecular mass of 151 kDa based on its primary sequence [START_REF] Claus | Molecular organization of selected prokaryotic S-layer proteins[END_REF]. Previous studies have observed a range of apparent molecular weights for SlaA from S. acidocaldarius DSM 639 up to 170 kDa [START_REF] Michel | The 2-D Crystalline Cell Wall of Sulfolobus acidocaldarius: Structure, Solubilization, and Reassembly[END_REF][START_REF] Grogan | Phenotypic characterization of the archaebacterial genus Sulfolobus: comparison of five wild-type strains[END_REF][START_REF] Selenska-Pobel | Magnetic Au nanoparticles on archaeal S-Layer ghosts as templates[END_REF]. Some early studies have even

P r e -P r i n t p r i o r t o p e e r r e v i e w

s e e f i n a l v e r s i o n d o i : 1 0 . 1 0 0 7 / s 0 0 7 9 2 -0 2 0 -0 1 1 7 9 -9

shown two high-molecular weight bands for the SlaA protein (140 kDa and 170 kDa) that were determined by their amino acid composition to be the same SlaA protein [START_REF] Michel | The 2-D Crystalline Cell Wall of Sulfolobus acidocaldarius: Structure, Solubilization, and Reassembly[END_REF]. This difference in the electrophoretic mobility of S. acidocaldarius S-layer proteins has been attributed to protein glycosylations, as well as their low solubility in acidic or neutral hot SDS [START_REF] Grogan | Organization and interactions of cell envelope proteins of the extreme thermoacidophile Sulfolobus acidocaldarius[END_REF]. In addition, S-layer proteins remain insoluble in the neutral pH Laemmli buffer containing SDS, and so prior to electrophoretic migration, S-layer proteins must be solubilized in a heated in alkaline buffer [START_REF] Grogan | Phenotypic characterization of the archaebacterial genus Sulfolobus: comparison of five wild-type strains[END_REF][START_REF] Selenska-Pobel | Magnetic Au nanoparticles on archaeal S-Layer ghosts as templates[END_REF][START_REF] Guan | N-glycosylation in the thermoacidophilic archaeon Sulfolobus acidocaldarius involves a short dolichol pyrophosphate carrier[END_REF]). S-layer proteins from other archaea have shown altered electrophoretic migration depending on culturing conditions including temperature, pH, and protein maturation events [START_REF] Konrad | Lipid modification of proteins in Archaea: attachment of a mevalonic acid-based lipid moiety to the surface-layer glycoprotein of Haloferax volcanii follows protein translocation[END_REF][START_REF] Claus | Molecular organization of selected prokaryotic S-layer proteins[END_REF]. While the reasons for these differences in electrophoretic mobility are still unknown, attention must therefore be paid when comparing results between studies using different strains of the same organism, as well as the culturing and protein solubilization conditions. This protocol allowed us to identify for the first time that the molecules aggregating with SlaB in strong association are GDGTs. The interaction between SlaB S-layer anchor protein and membrane tetraether lipids was shown to be highly resistant, even to heat and detergent treatments. Less SlaB is present than SlaA, even in whole cells, coherent with the current model for Sulfolobus S-layer structure with twice as much SlaA as SlaB [START_REF] Veith | Acidianus, Sulfolobus and Metallosphaera surface layers: structure, composition and gene expression[END_REF]).

There is likely some loss of SlaB during the initial steps of the extraction protocol due to strong attachment to membrane lipids. The final SlaB extract contained no other proteins detectable by SDS-PAGE or mass spectrometry, and no other lipids detectable using our HPLC-MS methods. Despite this, the SlaB proteins still eluted from the size exclusion chromatography column in the exclusion peak. This suggests that SlaB proteins strongly associate with one another, likely due to interactions between monomers, likely due to hydrophobic interactions between their N-terminal transmembrane domains. Very little is currently known about the role of SlaB. The SlaB anchor protein of S. solfataricus are predicted by its primary structure to be shorter than that of S. acidocaldarius, containing fewer beta sandwich domains [START_REF] Veith | Acidianus, Sulfolobus and Metallosphaera surface layers: structure, composition and gene expression[END_REF]. Correspondingly, the pseudoperiplasmic space created by the SlaB protein varies in height between the Sulfolobales, from 18 nm in Sulfolobus shibatae [START_REF] Baumeister | Principles of organization in eubacterial and archaebacterial surface proteins[END_REF] to 27 nm in S. acidocaldarius [START_REF] Kish | Preservation of Archaeal Surface Layer Structure During Mineralization[END_REF]. A recent study of ΔslaB mutants in Sulfolobus islandicus suggests that the S-layer anchoring role of SlaB can be partially complemented by another S-layer associated protein, designated M164_1049 [START_REF] Zhang | Cell Structure Changes in the Hyperthermophilic Crenarchaeon Sulfolobus islandicus Lacking the S-Layer[END_REF]. Another recent study in Sulfolobus solfataricus used CRISPR-mediated silencing of gene encoding SlaB to show that in addition to its anchoring role, SlaB is important in virus susceptibility, cell morphology, and cell division (Zink et al. 2019). These findings allow deeper study of the individual components of the essential cell envelope structure. This is of particular importance given the current absence of tertiary structures for archaeal S-layers, with the sole exception of the C-terminal domain of the Methanosarcina acetivorans S-layer protein [START_REF] Arbing | Structure of the surface layer of the methanogenic archaean Methanosarcina acetivorans[END_REF]. Extracted S-layers have been used in a range of biotechnology applications, from nanomaterials and bioremediation applications, to immobilized biocatalysts, biochips, and drug delivery [START_REF] Selenska-Pobel | Magnetic Au nanoparticles on archaeal S-Layer ghosts as templates[END_REF][START_REF] Bartolomé | Strong paramagnetism of gold nanoparticles deposited on a Sulfolobus acidocaldarius S layer[END_REF][START_REF] Sleytr | S-layers: principles and applications[END_REF]. Historically, S-layer protein purity is often verified in relation to other proteins by SDS-PAGE. Our results show the importance of analyzing for any attached lipids to ensure protein purity prior to any sensitive downstream applications. Taken together, the findings of this study will permit greater precision when isolating and studying

  Figures:

Figure 2 .Figure 4 .P

 24 Figure 2. Separation of SlaA and SlaB protein by gel filtration chromatography. After injection of the enriched S-layer protein sample obtained by hot SDS extraction (panel A) we obtain two peaks; the peak 1 is an exclusion peak containing SlaB, whereas the peak 2 corresponds to purified SlaA monomers. To make sure that no SlaA monomers were mixed with SlaB proteins in the excluded fraction, this fraction was then reinjected and separated by P r e -P r i n t p r i o r t o p e e r r e v i e w -s e e f i n a l v e r s i o n d o i : 1 0 . 1 0 0 7 / s 0 0 7 9 2 -0 2 0 -0 1 1 7 9 -9