Marien-Lorenzo Hanot
email: marien-lorenzo.hanot@umontpellier.fr.

NUMERICAL SOLUTION OF THE DIV-CURL PROBLEM BY FINITE ELEMENT EXTERIOR CALCULUS

Keywords: Finite Element, exterior calculus, Biot-Savart law, Hodge decomposition, Hodge-Dirac complex, de Rham complex, mixed element AMS subject classifications. 35F15, 65N30, 76B47, 76M10, 78M10

We are interested in the numerical reconstruction of a vector field with prescribed divergence and curl in a general domain of R 3 or R 2 , not necessarily contractible. To this aim, we introduce some basic concepts of finite element exterieur calculus [2] and rely heavily on recent results of P. Leopardi and A. Stern [9]. The goal of the paper is to take advantage of the links between usual vector calculus and exterior calculus and show the interest of the exterior calculus framework, without too much prior knowledge of the subject. We start by describing the method used for contractible domains and its implementation using the FEniCS library (see fenicsproject.org). We then address the problems encountered with non contractible domains and general boundary conditions and explain how to adapt the method to handle these cases. Finally we give some numerical results obtained with this method, in dimension 2 and 3.

1. Outline of the paper. In electromagnetism or fluid mechanics, one often encounter the problem of reconstructing a vector field with prescribed divergence and curl. As is well-known, the Biot-Savart law allows the reconstruction of a solenoidal vector field u in the whole space R d from its curl. u(x, t) = R d K(x, y) × curl u(y, t) dy with K(x, y) = 1 4π

x-y |x-y| 3

if d = 3 1 2π x-y |x-y| 2 if d = 2.
However it is a singular integral in the whole space and is not easily converted to a numerical algorithm, nor suitable in a bounded domain. The div-curl problem (2.1) has been addressed theoretically e.g. in [START_REF] Dacorogna | Existence and regularity of solutions of dω = f with Dirichlet boundary conditions[END_REF] or [START_REF] Auchmuty | L 2 -well-posedness of 3D div-curl boundary value problems[END_REF]. The numerical computation of the solution in a bounded domain is less documented, see e.g. [START_REF] Alonso Rodríguez | Simple finite element schemes for the solution of the curl-div system[END_REF] for a finite element solution using divergence free elements and an coercive variational form. The purpose of this paper is to show that the framework of differential forms and exterior calculus greatly simplifies both the theory and the finite element solution. In section 2 we show how the div-curl problem is related to the classical Helmholtz decomposition. In section 3 we introduce the main tools of exterior calculus. In section 4 we give a natural weak formulation of the div-curl problem, introduced in [START_REF] Leopardi | The Abstract Hodge-Dirac Operator and Its Stable Discretization[END_REF], which is well posed. In section 5 we detail the mixed finite elements compatible with the weak formulation. The implementation within the unified form language [START_REF] Logg | Automated solution of differential equations by the finite element method[END_REF] is sketched in section 6. The particular case of the 2D case is detailed in section 7. The case when the domain is not contractible is addressed in section 8. In section 9 we describe the cases of natural and essential boundary conditions. We close the paper with numerical tests in 2 and 3 dimensions in section 10.

Helmholtz decomposition.

Let Ω be a bounded domain in R d , the div-curl problem consists in finding a vector field with a prescribed divergence and curl. For g a scalar field and f a vector field we seek u such that (2.1)

∇ • u = g in Ω, ∇ × u = f in Ω.
Of course, one must add boundary conditions and specify the regularity of these fields. The existence of a solution is not guaranteed, indeed some vector fields are not written as the curl of other vector fields. The uniqueness is also an issue, indeed depending on the domain topology there may exist non trivial fields with vanishing divergence and curl, the so called harmonic fields. Theses problem make the elaboration of a stable scheme for the numerical solution quite complicated, because one must make sure that the fields f and g are compatible and in some way filter out the harmonic fields. Problem (2.1) is very close to the Helmholtz decomposition. The Helmholtz decomposition is central in vector calculus. In its classical formulation it states that any field of R 3 which is sufficiently smooth and decreases sufficiently fast at infinity can be decomposed into the sum of a gradient and a curl. The problem is then to compute this decomposition. For a given vector field F, find a vector potential A and a scalar potential φ such that (2.2)

F = -∇φ + ∇ × A.
Although this is not equivalent to the equations (2.1), we can see the relation by taking u = A in (2.2) and f = F. Thus, the problem of the absence of compatibility of f in (2.1) can be solved by looking simultaneously for a vector potential u and a scalar potential φ as described by (2.2). The general idea is then to couple these problems to obtain a well-posed system.

The classical Helmholtz decomposition (2.2) assumes that functions are smooth in the whole space. For bounded domains, we have the following standard result, see e.g. [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF] or chapter 9 of [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF]. Let Ω be a bounded, simply-connected, Lipschitz domain in R d , for any F ∈ (L 2 (Ω)) d there exists φ ∈ H 1 (Ω) and A ∈ H(curl, Ω) such that the decomposition (2.2) is valid. Of course, there is no uniqueness, since one can add constants to φ and any gradient field to A. Prescribing boundary conditions, we get the L 2 -orthogonal decompositions, proved e.g. in Arnold [START_REF] Arnold | Finite Element Exterior Calculus[END_REF]:

(L 2 (Ω)) d = ∇(H 1 (Ω)) ∇ × (H 0 (curl, Ω)) (2.3) (L 2 (Ω)) d = ∇(H 1 0 (Ω)) ∇ × (H(curl, Ω)). (2.4)
We recall the standard definitions:

H 1 0 (Ω) = {u ∈ L 2 (Ω); ∇u ∈ L 2 (Ω) d , u = 0 on ∂Ω}, H 0 (curl, Ω)) = {u ∈ L 2 (Ω) d ; curl(u) ∈ L 2 (Ω) d , u × n = 0 on ∂Ω}, n
being the unit outer normal to the boundary ∂Ω.

Exterior calculus.

Exterior calculus is based on differential forms, which are applications that at any point of the domain associate an alternating multilinear form. In particular the 0-forms are simply functions and the 1-forms can be seen as vector fields (by identifying linear forms with vectors, thanks to the usual inner product ). These spaces are endowed with a natural inner product (built from the one on vectors), which can be integrated on the whole domain. This allows us to define Hilbert spaces, in particular we will note L 2 Λ k (Ω) the set of L 2 -integrable k-forms. As alternating multilinear applications, it is possible to define an operation between a k-form α and an l-form β denoted by the wedge product ∧ giving a k + l-form α ∧ β = -β ∧ α. This operation is pointwise and allows to define a basis of the whole algebra from a basis of 1-forms. In dimension 3 a local basis is given in the table 3.1.

Space Λ 0 (Ω) Λ 1 (Ω) Λ 2 (Ω) Λ 3 (Ω) Basis 1
dx, dy, dz dy ∧ dz, -dx ∧ dz, dx ∧ dy dx ∧ dy ∧ dz Table 3.1: Basis of the exterior algebra on a domain of R 3 .

These forms are also equipped with the exterior derivative operator d. This operator acts globally on the exterior algebra d :

k=0 HΛ k (Ω) → k=0 HΛ k (Ω) but it is
often interesting to look at its action restricted to each space, where it will respect a notion of degree, specifically :

d : HΛ k (Ω) → HΛ k+1 (Ω) .
The operator d is defined as the usual differential on the 0-forms, by the property d • d = 0 and extended on the other forms degree by d(α ∧ β) = dα ∧ β + (-1) k α ∧ dβ for a k-form α. This differential operator is of course not defined on all L 2 -forms but only on a dense subset. As for the Sobolev spaces, we then consider the subset of L 2 -forms such that their exterior derivative is again L 2 . We note this set HΛ(Ω).

One can see a great similarity with vector calculus, for instance d • d = 0 corresponds to the well known identity curl(grad) = 0 and div(curl) = 0. This similarity is deep since there is a natural identification between vector fields and differential forms, the identification commutating with these differential operators. The identification is depicted by diagrams 3.2 for dimension 2 and 3.1 for dimension 3. We use the natural shortcut dx → x, dy → y, dz → z to mean that the one form a dx + b dy + c dz is identified with the vector a • e x + b • e y + c • e z = (a, b, c). In the same way, for 2-forms dy ∧ dz → x, dz ∧ dx → y, dx ∧ dy → z means that a dy ∧ dz + b dz ∧ dx + c dx ∧ dy is identified with the vector (a, b, c), whereas for 3-forms dx ∧ dy ∧ dz → 1 simply means that the 3 form a dx ∧ dy ∧ dz is identified with the scalar a. In dimension 2, the same conventions apply. We can notice that there are two possible identifications in dimension 2. We will refer to the first one in which ∇× appears as the curl identification and to the other one (in which ∇• appears) as the divergence identification. We refer again to [START_REF] Arnold | Finite Element Exterior Calculus[END_REF] for a thorough exposition of this so called "proxy" identifications.

Advantages of using exterior calculus are twofold. First a lot of work has already been done for the discretization of these spaces and operator d (see [START_REF] Arnold | Finite Element Exterior Calculus[END_REF], [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF]) as we will see in Section 5. Second is that the previously described vector -differential form identifications makes the nature of our operators clearer. Instead of having three different operators ∇, ∇• and ∇× we see that they become after identification a single operator d applied to different spaces. This allows us to unify the different types of "Helmholtz decomposition" (called Hodge decomposition in the context of exterior calculus) into a single one and allows us to understand how to find stable formulations of the problems (2.1) and (2.2).

Before stating this decomposition, let us introduce the adjoints of the operator d. Indeed, as we can see in the diagram 3.1, after identification the image of ∇ and that

H 1 (Ω) H(curl, Ω) H(div, Ω) L 2 (Ω) HΛ 0 (Ω) HΛ 1 (Ω) HΛ 2 (Ω) HΛ 3 (Ω) 1 x y z x y z 1 1 dx dy dz dy ∧ dz -dx ∧ dz dx ∧ dy dx ∧ dy ∧ dz ∇ ∇× ∇• d d d
Fig. 3.1: Identification between vectors and forms on a domain Ω ⊂ R 3 . of ∇× do not belong to the same space. This is perfectly normal because these two operators will not occur in the same decomposition, ∇ will occur in the decomposition of 1-forms and ∇× in that of 2-forms. The operators completing these decompositions are the adjoints of these operators. We define the codifferential operator δ as the adjoint operator of d. Since it is not defined on the whole L 2 Λ(Ω), we note its domain Ḣ Λ(Ω). Explicitly for Ω a bounded domain of R 3 , the adjoint of (∇, [START_REF] Arnold | Finite Element Exterior Calculus[END_REF] for details. Thus under the identification in fig. 3.1 we have

H 1 (Ω ) H(curl, Ω ) L 2 (Ω ) HΛ 0 (Ω ) HΛ 1 (Ω ) HΛ 2 (Ω ) 1 x y 1 1 dx dy dx ∧ dy ∇ ∇× d d H 1 (Ω ) H(div, Ω ) L 2 (Ω ) HΛ 0 (Ω ) HΛ 1 (Ω ) HΛ 2 (Ω ) 1 x y 1 1 dy -dx dx ∧ dy ∇ ⊥ ∇• d d
H 1 (Ω)) is (-∇•, H 0 (div, Ω)), the ad- joint of (∇×, H(curl, Ω)) is (∇×, H 0 (curl, Ω)) and the adjoint of (∇•, H(div, Ω)) is (-∇, H 1 0 (Ω)),
(3.1) (-∇•, H 0 (div, Ω)) = (δ, Ḣ Λ 1 (Ω)), (∇×, H 0 (curl, Ω)) = (δ, Ḣ Λ 2 (Ω)), (-∇, H 1 0 (Ω)) = (δ, Ḣ Λ 3 (Ω)).
Let Ω is a bounded, contractible, Lipschitz domain, to the two Helmholtz decompositions (2.3)-(2.4) correspond the two following Hodge decompositions, depending on the choice to identify (L 2 (Ω)) 3 with L 2 Λ 1 (Ω) or with L 2 Λ 2 (Ω). They are respec-tively given by

(3.2) L 2 Λ 1 (Ω) = d(HΛ 0 (Ω)) δ( Ḣ Λ 2 (Ω)), (3.3) L 2 Λ 2 (Ω) = δ( Ḣ Λ 3 (Ω)) d(HΛ 1 (Ω)).
Considering Λ 4 (Ω) ≡ {0} we can also extend the decomposition to the 3-forms:

L 2 Λ 3 (Ω) = d(HΛ 2 (Ω)).
This reflects the fact that

∇ • (H(div, Ω)) = L 2 (Ω).
Finally, if we try to extend the formula to 0-forms we will encounter a problem. Indeed, the constant functions are not in the range of (∇•, H 0 (div, Ω)). Moreover the constant functions are exactly the kernel of ∇ and are orthogonal to the range of (∇•, H 0 (div, Ω)), which is a generic fact that we will develop in Section 8 generalizing the algorithm to general domains. We will call such function harmonic (0-)forms and note their set H 0 ⊂ L 2 Λ 0 (Ω) (the set of constant 0-forms). The Hodge decomposition is then

L 2 Λ 0 (Ω) = δ( Ḣ Λ 1 (Ω)) H 0 .
We can gather these four decompositions into a decomposition of the total space

L 2 Λ(Ω) = 3 k=0 L 2 Λ k (Ω)
and on this domain, the Hodge decomposition is given by

(3.4) L 2 Λ(Ω) = d(HΛ(Ω)) δ( Ḣ Λ(Ω)) H 0 .
4. Weak formulation. The Hodge decomposition (3.4) gives the quasi invertibility of the operator d + δ : HΛ(Ω) ∩ Ḣ Λ(Ω) → L 2 Λ(Ω). We mean by that the following statement.

Theorem 4.1. Let f ∈ L 2 Λ(Ω), there exists w ∈ HΛ(Ω) ∩ Ḣ Λ(Ω) and h ∈ H 0 such that f = (d + δ)w + h.
Proof. Let f be in L 2 Λ(Ω), from the Hodge decomposition 3.4, there exist u ∈ HΛ(Ω), v ∈ Ḣ Λ(Ω) and h ∈ H 0 such that

f = du + δv + h.
Let us apply the Hodge decomposition to u and v. We can write u = dm + δp + q and v = dr + δs + t, where m, r belong to HΛ(Ω), p, s belong to Ḣ Λ(Ω) and q, t are in H 0 . Now take w = δp + dr. Using d 2 = δ 2 = 0 and the fact that q and t are harmonic forms, we can compute

(d + δ)w = dδp + δdr = d(dm + δp + q) + δ(dr + δs + t) = du + δv. So we have proved that f = (d + δ)w + h. Furthermore we have that δp ∈ L 2 Λ(Ω) and dm ∈ L 2 Λ(Ω), thus w ∈ L 2 Λ(Ω). Now dw = dδp = du belongs to L 2 Λ(Ω) and δw = δdr = δv belongs to L 2 Λ(Ω), hence w ∈ HΛ(Ω) ∩ Ḣ Λ(Ω).
This operator is called the Hodge-Dirac operator and is thoroughly studied in [START_REF] Leopardi | The Abstract Hodge-Dirac Operator and Its Stable Discretization[END_REF]. The non-invertibility of this operator comes from the space of harmonic forms H = H 0 , where we omit the superscript from now on, being both in the kernel of the operator and orthogonal to its range. This can be circumvented by using a space orthogonal to the harmonic forms.

Hence the primal formulation of (2.1) becomes:

Given f ∈ L 2 Λ(Ω) ∩ H ⊥ , find u ∈ HΛ(Ω) ∩ Ḣ Λ(Ω) ∩ H ⊥ such that ∀v ∈ HΛ(Ω), (4.1) du, v + δu, v = f, v .
By orthogonality arguments we can see that testing with v ⊥ H does not matter. Since d is adjoint to δ we can see that (4.1) is equivalent to:

Given f ∈ L 2 Λ(Ω) ∩ H ⊥ , find u ∈ HΛ(Ω) ∩ H ⊥ such that ∀v ∈ HΛ(Ω), (4.2) du, v + u, dv = f, v .
We remove the condition of orthogonality f ∈ H ⊥ by introducing a new pair of variables. The problem becomes:

Given f ∈ L 2 Λ(Ω), find u ∈ HΛ(Ω), p ∈ H such that ∀v ∈ HΛ(Ω), ∀q ∈ H, du, v + u, dv + p, v = f, v , (4.3) u, q = 0. (4.4)
We have added the component p, v to the equation (4.3), since by orthogonality we must have p = P H f (the orthogonal projection of f on H) and so we will effectively solve for (d + δ)u = f -P H f . Equation (4.4) ensures injectivity by imposing P H u = 0. Theorem 4.2. For any f ∈ L 2 Λ(Ω), there is a unique (u, p) ∈ HΛ(Ω) × H solution of (4.3)-(4.4). Moreover there exists c > 0 depending only on Ω such that u + du + p ≤ c f .

For the sake of completeness we give the proof, which follows closely the one given by Stern, [START_REF] Leopardi | The Abstract Hodge-Dirac Operator and Its Stable Discretization[END_REF]Theorem 6].

Proof. Let u ∈ HΛ(Ω) and p ∈ H, from the Hodge decomposition (3.4), there exist m ∈ HΛ(Ω), n ∈ Ḣ Λ(Ω) and o ∈ H such that u = dm + δn + o. Applying once again the Hodge decomposition to m, there exist r ∈ HΛ(Ω), s ∈ Ḣ Λ(Ω) and t ∈ H such that m = dr + δs + t. Take v = dδn + δs + p and q = o in (4.3)-(4.4). Noticing that dv = dδs = dm, du = dδn, using the Poincaré inequality δn ≤ c p dδn and using the orthogonality of the Hodge decomposition the equation reads:

(4.5) du, v + u, dv + p, v + u, q = dδn, dδn + u, dm + p, p + u, o = dδn 2 + dm 2 + p 2 + o ≥ 1 2 dδn 2 + 1 2c 2 p δn 2 + dm 2 + p 2 + o ≥ 1 2 du 2 + min(1, 1 2c 2 p ) u 2 + p 2 .
Moreover we can bound the HΛ-norm of (v, q) by the HΛ-norm of (u, p) with a Poincaré inequality:

(4.6)

v 2 + dv 2 + q 2 = dδn 2 + δs 2 + p 2 + dδs 2 + o 2 ≤ du 2 + c 2 p dm 2 + p 2 + dm 2 + o 2 ≤ (1 + c 2 p ) u 2 + du 2 + p 2 .
By the symmetry of the formulation this is enough to conclude with the Babuška-Lax-Milgram theorem.

Finally, let us translate this problem into the language of vector calculus in R 3 (the R 2 case being analogous with the appropriate definition of scalar and vector curl). First of all the unknowns are in fact a 4-tuple of fields, so we will write u = (u 0 , u 1 , u 2 , u 3 ), the subscript pertaining to the degree of the corresponding differential form. We will keep the notation H for the space of harmonic forms which is simply a vector space of dimension 1 containing the constant functions. The problem is then written:

Given (f 0 , f 1 , f 2 , f 3 ) ∈ L 2 (Ω) × (L 2 (Ω)) 3 × (L 2 (Ω)) 3 × L 2 (Ω), find u 0 ∈ H 1 (Ω), u 1 ∈ H(curl, Ω), u 2 ∈ H(div, Ω), u 3 ∈ L 2 (Ω), p ∈ H such that ∀v 0 ∈ H 1 (Ω), ∀v 1 ∈ H(curl, Ω), ∀v 2 ∈ H(div, Ω), ∀v 3 ∈ L 2 (Ω), ∀q ∈ H, (4.7) u 1 , ∇v 0 + p, v 0 = f 0 , v 0 , u 2 , ∇ × v 1 + ∇u 0 , v 1 = f 1 , v 1 , u 3 , ∇ • v 2 + ∇ × u 1 , v 2 = f 2 , v 2 , ∇ • u 2 , v 3 = f 3 , v 3 , u 0 , q = 0.
This weak formulation was introduced in [START_REF] Leopardi | The Abstract Hodge-Dirac Operator and Its Stable Discretization[END_REF], where its wellposedness is also established, thanks to the standard finite element exterior calculus tools [START_REF] Arnold | Finite Element Exterior Calculus[END_REF]. Natural boundary condition result from the weak formulation (4.7) hence a solution (u 0 , u 1 , u 2 , u 3 ) satisfies u 1 ∈ H 0 (div, Ω), u 2 ∈ H 0 (curl, Ω) and u 3 ∈ H 1 0 (Ω). Thus a solution of (4.7) satisfies (4.8)

-∇ • u 1 + p = f 0 , ∇u 0 + ∇ × u 2 = f 1 , ∇ × u 1 -∇u 3 = f 2 , ∇ • u 2 = f 3 .
with the natural boundary conditions (4.9)

u 1 • n = 0 on ∂Ω u 2 × n = 0 on ∂Ω u 3 = 0 on ∂Ω
Solving the problem (4.7) simultaneously solves two div-curl problems (2.1), computes two Helmholtz decompositions (of f 1 and f 2 ) and finds two functions with prescribed gradient. Assuming that the f i -functions are compatible (for the sake of simplicity only, otherwise the problem is solved for their orthogonal projection in the appropriate spaces) the problem div-curl (2.1) solved are:

∇ • u 1 = -f 0 ∇ × u 1 = f 2 u 1 • n = 0 on ∂Ω and ∇ • u 2 = f 3 ∇ × u 2 = f 1 u 2 × n = 0 on ∂Ω . Remark 4.3.
There are two differences between the two problems, first u 1 and u 2 do not follow the same boundary conditions, moreover as we shall see below in section 5, the discretization of the problem does not treat the differential and codifferential symmetrically. In particular, we have no error estimates for the convergence of the discrete codifferential ||δu -δ h u h || L 2 , see section 10.

5. Finite elements. The design for finite elements suitable for exterior calculus saw substantial progress recently, with the seminal work of [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF], [START_REF] Arnold | Finite Element Exterior Calculus[END_REF]. For the Hodge-Dirac problem (as well as for the Hodge-Laplacian problem), these elements can be realized in a very generic way as shown in the periodic table of the finite elements [START_REF] Arnold | Periodic table of the Finite Elements[END_REF].

The main properties of these elements are that they form a discrete subcomplex and admit bounded cochain projections. Being a discrete subcomplex means that the functions constructed on these elements belong to the domain of the exterior derivative (just as the functions are respectively included in H 1 , H(curl), H(div), L 2 ) and that their derivative (respectively their gradient, curl and div) are included in the functions of the next element. This last property allows to have properties exactly verified at the discrete level (for example discrete fields with exactly zero divergence). Bounded cochain projections are projections from continuous space to discrete space commuting with the exterior derivative. The boundedness for different norms ensures stability and accurate estimation of the error. These projections exist mainly as theoretical tools and their calculation is never performed in the numerical scheme.

Remark 5.1. Although the discrete exterior derivative operator d h is the same as the continuous d (more precisely its restriction on the discrete spaces), the discrete codifferential operator δ h has little to do with the continuous operator δ. Indeed, it is the adjoint of the same operator but on a different space. This is why we have removed any occurrence of δ from the formulation (4.1).

The space of harmonic forms H remains in this case the space of constant functions. Its determination can however be more complicated when using general domains or other boundary conditions. This problem is detailed in section 9.

Remark 5.2. The discrete spaces are then subspaces of continuous spaces and we use a conforming method. This will not always be the case for general domains where the space of discrete harmonic forms may be different from the continuous space, though they have the same dimensions, see [START_REF] Arnold | Finite Element Exterior Calculus[END_REF].

Although other types of meshes such as quadrilaterals are possible [START_REF] Arnold | Periodic table of the Finite Elements[END_REF], we focus on simplicial meshes. For each degree of forms, there are two families of piecewise polynomial elements indexed by their polynomial degree. Consider a simplicial mesh T (and denote the cells of the mesh by T ∈ T) and a polynomial degree r. The first family is the complete space of polynomials, differing by the type of continuity desired at the interfaces, it is given by (5.1)

P r Λ 0 (T) = {ω ∈ H 1 (Ω), ∀T ∈ T, ω |T ∈ P r (T, R)}, P r Λ 1 (T) = {ω ∈ H(curl, Ω), ∀T ∈ T, ω |T ∈ P r (T, R 3 )}, P r Λ 2 (T) = {ω ∈ H(div, Ω), ∀T ∈ T, ω |T ∈ P r (T, R 3 )}, P r Λ 3 (T) = {ω ∈ L 2 (Ω), ∀T ∈ T, ω |T ∈ P r (T, R)}.
The space P r (T, R k ) denotes the set of polynomials of degree r, defined on the domain T with value in R k . The second family are the so-called trimmed space with less degrees of freedom, it follows the same continuity conditions at the interfaces as the first one but uses only a subset of the polynomials, the exact definition of this set can be found in [START_REF] Arnold | Finite Element Exterior Calculus[END_REF], we will simply denote them here by P - r Λ k . This second family is then given by (5.2)

P - r Λ 0 (T) = {ω ∈ H 1 (Ω), ∀T ∈ T, ω |T ∈ P - r Λ 0 }, P - r Λ 1 (T) = {ω ∈ H(curl, Ω), ∀T ∈ T, ω |T ∈ P - r Λ 1 }, P - r Λ 2 (T) = {ω ∈ H(div, Ω), ∀T ∈ T, ω |T ∈ P - r Λ 2 }, P - r Λ 3 (T) = {ω ∈ L 2 (Ω), ∀T ∈ T, ω |T ∈ P - r Λ 3 }.
We always have P - r Λ 0 (T) = P r Λ 0 (T) and P - r Λ 3 (T) = P r-1 Λ 3 (T). The common names (as referred to in the table [START_REF] Arnold | Periodic table of the Finite Elements[END_REF]) for these elements in R 3 are

• Lagrange elements of degree r for P - r Λ 0 (T), • Nedelec's face elements of the first kind for P - r Λ 1 (T) and of the second kind for P r Λ 1 (T), • Nedelec's edge elements of the first kind for P - r Λ 2 (T) and of the second kind for P r Λ 2 (T),

• discontinuous Galerkin for P r Λ 3 (T). From these elements, we have to build a sequence of degree increasing forms. For each degree, we can take either the full polynomial element or the trimmed polynomial element. However, we must choose the appropriate polynomial degree, following a simple rule: if the next element is a complete polynomial we must go down one degree, if it is trimmed we keep the same polynomial degree. Thus for any degree r the sequences (5.3)

P - r Λ 0 (T) → P - r Λ 1 (T) → P - r Λ 2 (T) → P - r Λ 3 (T),
(5.4)

P r+3 Λ 0 (T) → P r+2 Λ 1 (T) → P r+1 Λ 2 (T) → P r Λ 3 (T)
are both correct.

The main difference between the two families comes from their approximation properties. We define the approximation error for a discrete space V h embedded in a continuous space V and a function u ∈ V by E(u) = inf v h ∈V h u-v h . Error estimates for schemes are often expressed in terms of these approximation errors. These errors depend on the size h of the cells in the mesh, and converge to 0 when h tends to 0. For these spaces we have in all cases an error estimate of the form E(u) ≤ C h l u H l with C a constant independent of h and l the order of convergence. For complete spaces of degree r the order of convergence of the approximation error (of a sufficiently regular function) is r + 1 and the order of convergence of its derivative is r. For trimmed spaces of degree r the order of convergence is r and the order of convergence of its derivative is also r.

Remark 5.3. To complete remark 4.3, we can see that in the discrete case another difference appears between the 1-forms and the 2-forms regarding the regularity of the solutions. According to whether we use Raviart-Thomas-Nedelec edge (resp. face) elements the discrete solution u h belongs to H(div, Ω) (resp. H(curl, Ω)). In the continuous case, the solution u belongs to both H(div, Ω) ∩ H(curl, Ω).

Let V 0 h → V 1 h → V 2 h → V 3
h be the chosen discrete space sequence and H h be the space of discrete harmonic forms, the discrete problem is then:

Given (f 0 , f 1 , f 2 , f 3 ) ∈ L 2 (Ω) × (L 2 (Ω)) 3 × (L 2 (Ω)) 3 × L 2 (Ω) find u h0 ∈ V 0 h , u h1 ∈ V 1 h , u h2 ∈ V 2 h , u h3 ∈ V 3 h , p h ∈ H h such that ∀v h0 ∈ V 0 h , ∀v h1 ∈ V 1 h , ∀v h2 ∈ V 2 h , ∀v h3 ∈ V 3 h , ∀q h ∈ H h , (5.5) u h1 , ∇v h0 + p h , v h0 = f 0 , v h0 , u h2 , ∇ × v h1 + ∇u h0 , v h1 = f 1 , v h1 , u h3 , ∇ • v h2 + ∇ × u h1 , v h2 = f 2 , v h2 , ∇ • u h2 , v h3 = f 3 , v h3 , u h0 , q h = 0.
The error estimates are calculated in [START_REF] Leopardi | The Abstract Hodge-Dirac Operator and Its Stable Discretization[END_REF]. Let K be the solution operator which takes f → u in (4.7) and π h be the bounded cochain projection mentioned at the beginning of section 5. We define η = (I -π h )K and µ = (I -π h )P H . In our case, both converge to 0, and in practice we can use their expressions computed in [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF], which for discrete spaces using polynomial degrees r give η = O(h) and µ = O(h r+1 ). The error estimate for u the continuous solution of (4.7) and u h the discrete solution of (5.5) gives:

(5.6)

d(u -u h ) E(du), u -u h E(u) + η(E(du) + E(p)), p -p h E(p) + µE(du).
Where the notation A B means that there is a constant C independent of u and h such that A ≤ CB. The equations (5.6) include all the components of u = (u 0 , u 1 , u 2 , u 3 ).

Implementation.

With these elements, the implementation of the problem (5.5) in unified form language (UFL), see [START_REF] Alnaes | Unified Form Language: A domain-specific language for weak formulations of partial differential equations[END_REF], is straightforward:

1 degree = 2 2 elemf0 = FiniteElement ( 'P ' , tetrahedron , degree ) 3 elemf1 = FiniteElement ( ' N1E ' , tetrahedron , degree ) 4 elemf2 = FiniteElement ( ' N1F ' , tetrahedron , degree ) 5 elemf3 = FiniteElement ( ' DG ' , tetrahedron , degree -1) 6 elemH = FiniteElement ( ' Real ' , tetrahedron , 0) 7 W = MixedElement ([ elemf0 , elemf1 , elemf2 , elemf3 , elemH ]) 8 9 ( u0 , u1 , u2 , u3 , uh ) = TrialF unctions ( W ) 10 ( v0 , v1 , v2 , v3 , vh ) = TestFunctions ( W ) 11 a1 = ( dot ( grad ( u0 ) , v1 ) + dot ( curl ( u1 ) , v2 ) + div ( u2 ) * v3 ) * dx 12 a2 = ( dot ( u1 , grad ( v0 ) ) + dot ( u2 , curl ( v1 ) ) + u3 * div ( v2 ) ) * dx 13 ah = uh * v0 * dx + u0 * vh * dx 14 a = a1 + a2 + ah Listing 1: Implementation of the variational formulation in ufl.

Our FEniCS codes and tests are available through the GitHub repository https:// github.com/mlhanot/divcurl solver.

Being a mixed finite element scheme, the assembled linear system is semidefinite, which for very large number of degrees of freedom could be an issue. However for contractible domains we obtained reliable results on a standard computer (16 GB RAM) with direct solvers up to several millions of degrees of freedom. For non contractible domains, the main workload is generated by the computations of a linear basis of the harmonic forms. Indeed this amounts to computing a basis of the null space of a matrix, which is trickier than solving a regular linear system.

7.

Problem in two dimensions. The exterior calculus framework is unified w.r.t. dimensions, hence the formulation (4.3)-(4.4) remains perfectly valid without any change (except that one needs Ω to be domain of R 2 instead of R 3 ).

The expression in the vector calculus formalism is however quite different, in particular there are two possibilities of identification. The two corresponding problems are then (respectively for the curl identification and the div identification):

Given (f 0 , f 1 , f 2 ) ∈ L 2 (Ω) × (L 2 (Ω)) 2 × L 2 (Ω), find u 0 ∈ H 1 (Ω), u 1 ∈ H(curl, Ω), u 2 ∈ L 2 (Ω), p ∈ H such that ∀v 0 ∈ H 1 (Ω), ∀v 1 ∈ H(curl, Ω), ∀v 2 ∈ L 2 (Ω), ∀q ∈ H, (7.1)
u 1 , ∇v 0 + p, v 0 = f 0 , v 0 , u 2 , ∇ × v 1 + ∇u 0 , v 1 = f 1 , v 1 , ∇ × u 1 , v 2 = f 2 , v 2 , u 0 , q = 0. Given (f 0 , f 1 , f 2 ) ∈ L 2 (Ω) × (L 2 (Ω)) 2 × L 2 (Ω), find u 0 ∈ H 1 (Ω), u 1 ∈ H(div, Ω), u 2 ∈ L 2 (Ω), p ∈ H such that ∀v 0 ∈ H 1 (Ω), ∀v 1 ∈ H(div, Ω), ∀v 2 ∈ L 2 (Ω), ∀q ∈ H, (7.2) u 1 , ∇ ⊥ v 0 + p, v 0 = f 0 , v 0 , u 2 , ∇ • v 1 + ∇ ⊥ u 0 , v 1 = f 1 , v 1 , ∇ • u 1 , v 2 = f 2 , v 2 , u 0 , q = 0.
The difference between the two is then the same as in the remark 4.3, namely the boundary conditions and the preferred operator.

The choice of finite elements must also be appropriate, although similar to the 3-dimensional case, one must be careful about the vector element chosen. The appropriate elements are (using the same semantics as for (5.1) and (5.2)):

• Lagrange elements for P - r Λ 0 (T), • Raviart-Thomas face (or edge) elements for P - r Λ 1 (T), • Brezzi-Doublas-Marini face (or edge) elements for P r Λ 1 (T),

• discontinuous Galerkin for P r Λ 2 (T). We must use edge elements for the problem (7.1) and face elements for the problem (7.2).

Remark 7.1. The problem (7.2) is in fact simply the problem (7.1) to which we apply a quarter rotation on the vector space. This means that we solve the problem in the space of differential forms and only then apply the identification of Fig. 3.2 to come back to the vector fields.

8. Non contractible domain and harmonic forms. When the domain is no longer contractible, the Hodge decomposition as given in (3.4) is no longer valid. As described in section 3 the problem comes from the appearing of harmonic forms, i.e. elements in the kernel of the Hodge-Dirac operator (equivalently in the kernel of the Hodge-Laplacian or fields f such that ∇ • f = 0 and ∇ × f = 0). This is not really an issue because we already had harmonic forms in (4.1) and we can treat them in the same way. However this poses the problem of determining the space of harmonic forms The theoretical aspect of the problem is solved thanks to the famous theorem of De Rham giving the isomorphism between the harmonic forms and the cohomology of the cochain complex. This cohomology has a very strong geometrical interpretation, its dimension is given by the Betti numbers. Thus in dimension 2 the number of harmonic 0-forms is the number of connected components of the domain, the number of harmonic 1-forms is the number of holes of the domain and there are no harmonic 2-forms. In dimension 3 the number of harmonic 0-forms is still the number of connected components, the number of harmonic 1-forms corresponds to the number of tunnels and the number of harmonic 2-forms corresponds to the number of vacuum bubbles, there is no harmonic 3-forms. This is illustrated in figure 8.1 showing the two harmonic 1-forms on a disk with two holes, and figures 8.2 and 8.3 showing the two harmonic 1-forms on a hollow torus.

Another important theorem gives the isomorphism between discrete and continuous harmonic forms (see [START_REF] Arnold | Finite Element Exterior Calculus[END_REF], chapter 5 and 7.6), so the dimension of the space of harmonic forms does not depend on the discretization or the elements chosen.

However the actual search for these forms is much more complicated. The isomorphisms only give their dimension number, and a heuristic idea of their shapes. To compute them, we will start from their definition as kernel of the Hodge-Dirac operator, which corresponds exactly to the kernel of the assembled matrix of the system (see [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF]). The problem of determining a basis of harmonic forms becomes in practice a problem of finding kernels of matrices. The dimension of these kernels is a useful information for many algorithms (for example to search for the smallest eigenvalues) and the idea of the solution shapes can give a good initial guess. However, the problem remains intricate because of the size of the linear systems.

In our numerical computations with FEniCS, we have achieved the best results with the numerical algebra library SLEPc, see https://slepc.upv.es/.

Once a basis of harmonic forms is determined, it is enough to add them to the space H in the scheme. The proofs are done in this general framework and still give the right estimates. 9. Boundary conditions. So far, we have not imposed any essential (Dirichlet) conditions on our spaces, so natural conditions have emerged. Although these conditions are sufficient to ensure wellposedness, they fix the degree of forms in which we look for our solution and we have noticed in the remarks 4.3 and 5.3 that these choices have an impact on the convergence and the regularity of our solutions. Specifically, natural boundary conditions u × n = 0 (resp. u • n = 0 ) imposes u to be a 2-form, i.e. u ∈ H(div, Ω) (resp. u to be a 1-form, i.e. H(curl, Ω)). If we wish that u h belongs to the space which does not correspond to the natural condition obtained, a simple way is to apply homogeneous Dirichlet conditions to all spaces. The sequence then becomes (9.1)

H 1 0 (Ω) H 0 (curl, Ω) H 0 (div, Ω) L 2 (Ω). ∇ ∇× ∇•
From a theoretical point of view, enforcing Dirichlet conditions on all spaces does not pose any problem, we just end up with another complex, dual of the first one and all the theorems still work with one difference: the number (dimension) of harmonic forms are inverted, so under these conditions there are no 0 harmonic forms and as many 3 harmonic forms (in dimension 3) as there are connected components. One must adjust the H space in the formulation, so replace p, v 0 by p, v 3 in (4.7).

The situation becomes more complicated when mixed conditions are applied (natural on some faces, essential on others). New harmonic forms can then appear even for simple domains as illustrated in figure 9.1.

10. Numerical application. To conclude we present some results computed with the weak form (5.5). In dimension 2 we take for reference function (10.1) on the unit square (0, 1) 2 , tailored to accommodate the boundary conditions, to have a non-trivial divergence and curl without being symmetric. In dimension 3 we take for reference function (10.2) on the unit cube (0, 1) 3 . The test with non-trivial harmonic forms is delicate, indeed the reference function will not be orthogonal to these forms and these forms are not known for most of the meshes. Thus we computed the rate of convergence in dimension 2 only, using periodic conditions on the edges (x = 0) and (x = 1) of the unit square with vanishing tangential components on the other edges because the harmonic forms are explicitly known in this case, see Fig. 9 We report the rate of convergence of the discrete solution towards the analytical solution, the rate of convergence of its differential (divergence or curl as the case may be) and the rate of convergence of its codifferential. The latter is calculated in two steps, first we make the orthogonal projection P of the function on an appropriate space (if the function has been defined on the face elements we project it on the corresponding edge elements and inversely). Then we compute its codifferential. Although there is no theoretical result on this convergence to our knowledge, we can observe in the last columns of the tables a convergence rate for the codifferential but with a loss of one degree with respect to the approximation property of the space of u h . Of course as reported in tables 10.1, 10.5, 10.7 and 10.9, when degree one polynomials are used no convergence of the codifferentials can be expected.

Using a contractible domain and the sequence of finite elements (5.3), we can see the convergence rates in dimension 2 for the divergence identification for polynomials of degree 1 and 2 in tables 10.1 and 10.2, for the curl identification in tables 10.5 and 10.6. In dimension 3 we can see the convergence rates for 1-forms in tables 10.7 and 10.8, for 2-forms in tables 10.9 and 10.10. We can see the convergence rates with the sequence of elements (5.4) of degree r = 0 and r = 1 in dimension 2 with the divergence identification in tables 10.3 and 10.4. Finally, we can see the rate of convergence when harmonic 1-forms are present in dimension 2 with the divergence identification using the sequence (5.4) of polynomial degree 2 in the table 10.11. We notice in tables 10.3 and 10.4 that using different polynomials degrees like in equation (5.4) we have no loss of convergence for the codifferential. We have no theoretical proof of this fact.

Fig. 3 . 2 :

 32 Fig. 3.2: Two possibles identification between vectors and forms on a domain Ω ⊂ R 2 .

Fig. 8 . 1 :

 81 Fig. 8.1: Harmonic 1-forms on a surface with two holes.

Fig. 8 . 2 :

 82 Fig. 8.2: First harmonic 1-form on a 3D hollow torus sliced for visualisation.

Fig. 8 . 3 :

 83 Fig. 8.3: Second harmonic 1-form on a 3D hollow torus sliced for visualisation.

Fig. 9 . 1 :

 91 Fig. 9.1: Harmonic 1-form on a contractible domain with mixed boundary conditions.

Table 10 .

 10 .1. 1: Convergence rates for the sequence (5.3) of degree 1 in 2-dimension with the divergence identification.

	(10.1)	u =	sin(3πx) cos(πy) sin(πy) cos(2πx)

Acknowledgements. The first author would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme GCS-Geometry, compatibility and structure preservation in computational differential equations, 2019 where work on this paper was initiated. This work was supported by EPSRC grant no EP/R014604/1.

References.


Delphine Herrera

Thibaud Chevalier

Didier Frot

Loïc Barré

Audrey Drelich

Isabelle Pezron

Christine Dalmazzone
email: christine.dalmazzone@ifpen.fr

Monitoring the formation kinetics of a bicontinuous microemulsion

Keywords: Microemulsions, Bicontinuous, Formation, Quantification, Structure, DLS, NMR, Micro-CT, SAXS

Hypotheses

The performance of bicontinuous microemulsions is usually assessed on the characteristics of the middle phase at equilibrium. However, applied to Enhanced Oil Recovery, such an evaluation would not be representative of the structure and composition of fluids in reservoir rocks. Studies on the properties of non-equilibrated microemulsions are still needed to better understand the formation of such complex systems, in particular to optimize input parameters of process simulation tools.

Experiments

For this purpose, we monitored the formation of a microemulsion from contact with the oil to equilibrium when no mixing or convection is provided. Non-destructive methods such as Nuclear Magnetic Resonance, Micro-Computed Tomography, Dynamic Light Scattering and Small Angle Xray scattering were used to extract the compositions, phase thicknesses, dynamics and structures of the system over time.

Finding

We found that the system gets structured into several layers over time that include the transient presence of an oriented semi-crystalline phase. The growth of the bicontinuous middle phase results from a progressive reorganization of the liquid crystal. The compositional and structural gradients, observed along the sample height, are correlated and linked to the corresponding structures of the phase diagram of the quaternary system. Equilibrium is reached after the total transfer of the liquid crystal into the bicontinuous phase.

Introduction

Since their discovery in 1959 [START_REF] Schulman | Mechanism of Formation and Structure of Micro Emulsions by Electron Microscopy[END_REF], microemulsions have been intensively studied due to their remarkable properties such as low viscosity, isotropic and thermodynamically stable phases as well as ultralow interfacial tensions between oil and water [START_REF]Handbook of microemulsion science and technology[END_REF] [START_REF] Danielsson | The definition of microemulsion[END_REF]. They are used today in a wide variety of applications in all fields of chemistry such as cosmetic, chemical, pharmaceutical, food and oil industries mainly for their high oil solubilization performance. Particularly in Enhanced Oil Recovery (EOR), the use of microemulsions appears to be a good alternative to extract the crude that is trapped by capillary forces into the porous reservoir rock. When injecting an optimal surfactant-rich formulation, a microemulsion should form after contact with the crude. When anionic surfactants are used, a co-surfactant (alcohol type) is often added to the formulation to enhance the interfacial film fluidity and allow the formation of the microemulsion. In most studies, the optimal formulation is determined according to the performances of the microemulsion at equilibrium. However, this evaluation may not be representative of the behavior of fluids in situ in the reservoir [START_REF] Unsal | Pore Scale Dynamics of Microemulsion Formation[END_REF]. In a real porous rock, two situations may occur: either the pores are interconnected, or they are dead-end pores. In the first case, convection processes ensure a natural mixing of fluids in situ leading to a quasi-instantaneous formation of the microemulsion. In the second case, no source of mixing is provided. Recently, Broens et al. [START_REF] Broens | Emulsification kinetics during quasi-miscible flow in dead-end pores[END_REF] studied such microemulsification in a microfluidics set up by modeling the dead ends pores of a reservoir. In this configuration, they confirmed that the formation of the microemulsion is only controlled by diffusion phenomena that take place in the pore space. By monitoring the penetration front of the microemulsion under fluorescence microscopy, they indicated that the kinetics of micro-emulsification obeys Fick's law. However, if matter exchanges are limited by diffusion, one can wonder if the bicontinuous microemulsion is properly formed in situ. In 2019, Unsal et al. [START_REF] Unsal | Imaging of compositional gradients during in situ emulsification using X-ray micro-tomography[END_REF] investigated the micro-emulsification of a bicontinuous microemulsion under flow in a model porous rock using quantitative X-ray microtomography. They directly proved the presence of local gradients of composition in the vicinity of the dead-end pores in the porous rock. However, current methodologies and analytical instruments are limited in accessing local properties and structure of fluids in these complex regions. To some extent, these data would be particularly useful for optimizing oil recovery processes and better understanding plugging issues during flow. We thus propose to work on the mechanisms of formation of a bicontinuous microemulsion, in bulk and under static conditions, when only diffusion processes are present. The challenge is to have access to local and quantitative information over time without affecting the diffusion phenomena occurring in situ in the system. In particular, analytical methods that allow quantitative kinetic studies of the formation of microemulsions are scarce. The usual chemical techniques such as Karl-Fischer or Hyamine titrations [START_REF] Fukumoto | Investigation on Physical Properties and Morphologies of Microemulsions formed with Sodium dodecylbenzenesulfonate, Isobutanol, Brine, and Decane, Using Several Experimental Techniques[END_REF][START_REF]Micellar Solutions and Microemulsions: Structure, Dynamics, and Statistical Thermodynamics[END_REF] as well as gas chromatography [START_REF] Huie | Partition of n-Butanol among phases and solubilization ability of Winsor type III Microemulsions[END_REF][START_REF] Huie | Solubilization Behavior of Organic Mixtures in Optimum Winsor Type III Microemulsion Systems of Sodium dodecylsulfate[END_REF] are destructive and relatively time-consuming for kinetic perspectives. In addition, these methods do not allow simultaneous and local quantification of a multi-component system. In this objective, we benefited from the complementarity of Low-Field Nuclear Magnetic Resonance (NMR), X-Ray Micro-Computed Tomography (Micro-CT), Dynamic Light Scattering (DLS) and Small Angle X-Ray Scattering (SAXS) to monitor the formation of a bicontinuous microemulsion from phases contact to equilibrium. The studied microemulsion is a quaternary system containing toluene/NaCl brine/butan-1-ol/Sodium dodecylsulfate [START_REF] Bellocq | Three-Dimensional Phase Diagram of the Brine-Toluene-Butanol-Sodium dodecylsulfate System[END_REF] at optimal formulation. The in-house NMR method was developed to obtain composition profiles of oil, brine and co-surfactant along the sample height [START_REF] Herrera | Quantification of microemulsion systems using low-field T1-weighted imaging[END_REF]. Micro-CT enabled precise acquisitions of phase thicknesses and volume fraction of brine in the sample. Finally, we correlated the compositional information collected from these techniques with SAXS and DLS measurements which provided respectively local and dynamic information on the structure. The methods implemented in this work are non-destructive and the measurements are fast enough to allow kinetic monitoring.

Experiments

Products

The surfactant, Sodium dodecylsulfate, (dry weight > 98.5%, water <1.5%) was supplied by Alfa Aesar and was used without further purification. Toluene was from VWR (AnalaR Normapur), butan-1-ol and sodium chloride NaCl were from Merck (grade analysis). Milli-Q water was used (HX7000 from Merck, Resistivity 10-15 MΩ.cm; Conductivity < 0.2 µS/cm at 25°C).

Microemulsions systems

The Winsor III system studied is made of NaCl brine, Sodium dodecylsulfate "SDS" (surfactant), butan-1-ol (co-surfactant) as aqueous phase and toluene as oil phase. Compositions were calculated to satisfy a water-to-oil volume ratio WOR = 1; co-surfactant = 5.5 % v/v. of the total system; cosurfactant/Surfactant weight ratio = 2. More details on the formulation steps can be found in Supplementary Material Part I. Experiments were carried out on a microemulsion at optimal formulation (Hydrophilic Lipophilic Deviation = 0) [START_REF] Salager | Partitioning of Ethoxylated Octylphenol Surfactants in Microemulsion-Oil-Water Systems: Influence of Temperature and Relation between Partitioning Coefficient and Physicochemical Formulation[END_REF] corresponding to an aqueous phase salinity of 54 g/L of NaCl. The microemulsion was directly formulated in 1.5 mm-diameter quartz capillaries for SAXS (50 µL) and in 1.8 cm diameter glass tubes for NMR (3 mL), Micro-CT and DLS (10 mL).

Kinetics studies

Kinetic studies were performed by NMR, Micro-CT and DLS. SAXS was used punctually to characterize the structures encountered during the process. We first independently formulated the aqueous phase of the microemulsion as described above. Toluene was added at a controlled speed of 126.8 mL/h via an automatic syringe pump on the aqueous phase without any stirring or mixing. After Toluene addition, the vials were sealed with caps and Teflon rubber to avoid evaporation. During NMR experiments, the samples were inserted into the probe and the temperature was set at 30°C by temperature-controlled air stream. While in the DLS, Micro-CT and SAXS tests, the experiments were carried out at room temperature (20°C) as no temperature control device is present on these experimental setups. To verify if any change of structure occurs at a warmer temperature (between 20 and 30°C), the equilibrated microemulsion was left to rest 1 day in an oven at 40°C. We did not notice any change in the thicknesses of the triphasic system nor any change in the aspect of the intermediate phase. Thus we assume that the temperature difference between the setups environment does not significantly affect the structure and composition of each phase in the sample which are still comparable. The measurements along the tube height were performed under completely static conditions until equilibrium was reached. Once equilibrium is reached, the system is triphasic with a bicontinuous microemulsion in the middle phase and excess oil and aqueous phases on both sides as illustrated in Figure 1 [14] [START_REF] Duvail | Mesoscopic modelling of frustration in microemulsions[END_REF]. We define the equilibrium state as the point where no more changes in the interfaces positions are noticed and constant properties are found in each phase with the characterization techniques used. Low-Field NMR acquisitions were carried out on a 20 MHz benchtop MQR spectrometer from Oxford instruments equipped with a probe of 1.8 cm outer-diameter. The set up allows to vertically scan a tube over 2 cm in height corresponding to a 3 mL sample. The temperature was controlled at 30°C during the experiment via the injection of compressed air into the probe. For one measurement, the acquisition time is 1 min 40 s. The kinetics monitoring was stopped after 7 days of stabilization. We developed [START_REF] Herrera | Quantification of microemulsion systems using low-field T1-weighted imaging[END_REF] a methodology that allows the simultaneous quantification of toluene, brine and butan-1-ol along the tube height. It is worth mentioning that the same method could potentially be applied to any multicomponent, dark/opaque and viscous system insofar as the sample is liquid and contrasts are present between compounds. The accuracy of the method does not allow the quantification of small quantities lower than 4 % v/v in the sample. To allow data treatment, a hypothesis had to be made on SDS partition in the tube: we assumed the surfactant is entirely and homogeneously distributed in the intermediate phase(s) between excess oil and brine phases at equilibrium. This hypothesis has been verified in the literature on the same systems by carrying out Hyamine titrations or atomic adsorption tests in the excess aqueous phase and middle phase [START_REF] Fukumoto | Investigation on Physical Properties and Morphologies of Microemulsions formed with Sodium dodecylbenzenesulfonate, Isobutanol, Brine, and Decane, Using Several Experimental Techniques[END_REF][START_REF] Bellocq | Three-Dimensional Phase Diagram of the Brine-Toluene-Butanol-Sodium dodecylsulfate System[END_REF]. It is shown that 0 to 1 % w/w of surfactant remains in the excess water phase thus the surfactant is quasi-totally distributed in the middle phase. When no intermediate phase can be identified, we assumed that SDS was distributed throughout the aqueous phase. To quantify the other three compounds of the system the principle is based on 1D-imaging using weighted T1relaxometry. Detailed information on the calculations of volume fractions from T1-relaxation profiles can be found in Supplementary Material Part II. The composition profiles resolution along the height of the tube is 100 µm.

Phase behavior and quantification by X-ray Micro-Computed Tomography

The X-Ray Micro-Computed Tomography (Micro-CT) tests were conducted on an EasyTom 150-160 from RX Solutions in 10 mL glass tubes of 1.8 cm-diameter at ambient temperature (~20°C). The principle is based on the transmission of X-rays radiation through the sample. In the literature, the technique is commonly used to quantify the brine from oily compounds (here, butan-1-ol and toluene) because of a strong contrast in the attenuation coefficients [START_REF] Unsal | Imaging of compositional gradients during in situ emulsification using X-ray micro-tomography[END_REF][START_REF] Scheffer | Enhancement of oil recovery by emulsion injection: A pore scale analysis from X-ray micro-tomography measurements[END_REF][START_REF] Withers | Quantitative X-ray tomography[END_REF]. Practically, radiographies of the sample were acquired over time with a duration of 5 s/radiography as displayed in level of greys in Figure 5. The raw transmitted intensities were collected and corrected with an empty tube signal as refernce. The calculation of the brine volume fractions along the tube is deduced from Equation 1. However, no significant contrast was observed between the references [5 % w/w SDS in water] and [water]. Therefore, the presence of SDS solubilized in brine within the different phases of the microemulsion could not be detected by this technique. Moreover, the Micro-CT analysis provides precise information on the phase thicknesses and interfaces positions with a partial pixel resolution of 72 µm.

With 𝑥 (ℎ): the volume fraction of brine at height h in the sample, 𝐼 and 𝐼 : the respective transmitted intensities of separate brine reference and toluene reference (or butan-1-ol reference as the same transmitted intensities were measured for toluene and butan-1-ol).

Dynamic Light Scattering (DLS)

Measurements were performed using a DLS apparatus, Vasco-flex from Cordouan Technologie. The setup is equipped with a laser diode as light source (λ = 656 nm) and a single mode fiber connected to a photon counting detector. Thanks to the flexibility of the fiber, the sample was directly analyzed in its container (10 mL and 1.8 cm-diameter glass tube) so that the internal structure of the microemulsion was preserved. The diffusion measurements were conducted at an angle 𝜃 of 90° at room temperature (20°C). Data were collected and analyzed by the NanoQ software. The apparatus measures fluctuations of intensity versus time and plots the autocorrelation function G(t) of the scattered intensity. For random processes, the autocorrelation function is properly fitted by an exponential decay in the form written in Equation 2If the dispersed objects are considered as spherical hard spheres having a Brownian motion,

𝑥 (ℎ) = 𝑙𝑛(𝐼 ) -𝑙𝑛 𝐼 (ℎ) 𝑙𝑛(𝐼 ) -𝑙𝑛(𝐼 ) Equation 1
the Stokes-Einstein relation (Equation 3) is commonly used to determine the particles size from the measurement of the characteristic relaxation time 𝜏. In the case of a bicontinuous system or non-Brownian motion of particles, the Stokes-Einstein relation is no longer applicable and the interpretation of the relaxation times becomes more complex. In the literature, DLS studies on bicontinuous structures are scarce but one can interpret DLS relaxation times of bicontinuous microemulsions as the ability to create/break connections between oil and water domains that is to say 𝜏 is probably relative to the rate of matter exchanges through the interfacial film of oil and water domains [START_REF] Fukumoto | Investigation on Physical Properties and Morphologies of Microemulsions formed with Sodium dodecylbenzenesulfonate, Isobutanol, Brine, and Decane, Using Several Experimental Techniques[END_REF][START_REF] Sharma | Incorporation of Melittin Enhances Interfacial Fluidity of Bicontinuous Microemulsions[END_REF]. Low relaxation times could be potentially linked to highly interconnected interfaces. We thus noted Dapp, the apparent diffusion coefficient that stems from Equation 3 when referring to the bicontinuous structures or a non-Brownian motion of particles.

𝐺(𝑡) = 𝐴 𝑒𝑥𝑝(-2𝐷𝑞 𝑡) and 𝜏 = ² Equation 2
With 𝐴 a normalization constant, 𝜏 the characteristic relaxation time of the process in µs, 𝑞 = 𝑠𝑖𝑛 (𝜃 2 ⁄ ) the wave vector in m -1 , and D the diffusion coefficient in m².s -1 .

𝐷 = 𝑘𝑇 3𝜋𝐷 𝜂 Equation 3
With k the Boltzman constant in kg.m 2 .s -2 .K -1 , T the temperature in K, 𝜂 the dynamic viscosity of the continuous medium in Pa.s, 𝐷 the hydrodynamic diameter of the dispersed object in m.

Small Angle X-Ray Scattering (SAXS)

SAXS measurements were conducted on a Xeuss 2.0 set up from Xenocs on the Winsor III system. The microemulsions were directly formulated into quartz capillaries (diameter of 1.5 mm) to avoid sampling. We first inserted 25 µL of a pre-formulated aqueous phase using a Hamilton syringe. 23 µL of toluene was added via an iron rod above the aqueous phase to prevent any mixing in the capillary. The capillary was sealed by a drop of silicon paste. We studied two samples that were formulated in two different capillaries: a non-equilibrated system obtained 3 days after toluene addition and the equilibrated bicontinuous microemulsion measured at 7 days after toluene addition. It is worth noting that the addition speed of Toluene in capillaries cannot be controlled and implies that the equilibration times are not comparable with DLS, NMR or Micro-CT experiments even if the same processes are observed. The sample-to-detector distance was 600 mm allowing a q-range from 0.01 to 0.40 Å -1 . The X-ray beam (λ=1.54 Å) was collimated by two slits: 0.7 mm² and 0.4 mm². The calibration of the q-range was performed using standard silver behenate AgBe as reference. The intensity calibration was carried out with Lupolen only for the equilibrated sample. The cross-section intensities are given in cm -1 for this sample. Practically, the transmitted intensities were scanned from bottom to top in order to locate the intermediate phases. Accordingly, scattering measurements were performed at various heights. For the non-equilibrated sample, we confirmed that the system did not evolve during the experiment by carrying out another transmission scan after the experiment that overlaps the previous scan.

In the case of the bicontinuous microemulsion system, scattering was found to be isotropic and the azimuthally averaged curves I(q) were fitted using the well-known Teubner-Strey model [START_REF] Teubner | Origin of the scattering peak in microemulsions[END_REF]20] widely used for microemulsions SAXS data interpretation in medium q range (more information on the model can be found in Supplementary Material Part III). It gives access to the periodicity of the structure 𝑑 and the correlation length 𝜉 (i.e. the characteristic decay length of periodic order). On the SAXS curve, the position qmax of the characteristic peak of the microemulsions is related to d = 2π/qmax and its width is related to ξ: the longer the correlation length is, the thinner the peak will be [START_REF] Oberdisse | Structure, interfacial film properties, and thermal fluctuations of microemulsions as seen by scattering experiments[END_REF]. From these two parameters the amphiphilicity factor fa can be calculated (Equation 4). The amphiphilicity factor fa indicates a scale of the extent of ordering of the molecules at the interfacial film within Winsor III systems. It goes from +1 with completely disordered water/oil interfaces to -1 with well-ordered interfaces leading to lamellar phase.

𝜉 𝑑 = 1 4𝜋² 1 -𝑓𝑎 1 + 𝑓𝑎

Equation

In the case of an oriented phase [START_REF] Davidson | The measurement of the Nematic Order Parameter by X-ray Scattering Reconsidered[END_REF][START_REF] Maier | Eine einfache molekular-statistische Theorie der nematischen kristallinflüssigen Phase[END_REF][START_REF] Gr | Why is the Maier-Saupe theory of nematic liquid crystals so successful?[END_REF], we adjusted the azimuthal profiles by a Maier-Saupe function (Equation 5) to have access to the Hermann or order parameter S (Equation 7) corresponding to the second moment of the normalized Maier-Saupe function (Equation 6)When S is equal to 1, the structure is perfectly ordered, when it is equal to 0 the structure is random. The fit also provides the angle of orientation α of the scattering objects.

𝐼(𝜓) = 𝐴 + 𝐵exp(𝑚cos²(𝜓 -𝛼)) Equation 𝑓 norm (𝛽) = 𝐴 + 𝐵exp(𝑚cos²(𝛽)) 2𝜋 ∫ (𝐴 + 𝐵exp(𝑚cos²(𝛽)))sin(𝛽)𝑑𝛽 Equation 𝑆 = 3cos²( -𝛽) -1 2 sin = 𝑑𝜙 𝑓 norm (𝛽) 3cos²(𝛽) -1 2 sin(𝛽)𝑑𝛽

Equation

Where A and B are normalization constants; m is the orientation parameter accounting for the liquid crystal interactions; 𝜓 is the azimuthal angle; α the orientation angle of the liquid crystals; β = 𝜓α.

Results and discussion

Equilibrium state characterization

The compositions of the final system after equilibration calculated from Micro-CT and NMR profiles are given in Table 1. A good agreement is found between the two techniques showing that the kinetics monitoring could be carried out. Micro-CT is more reliable than NMR in assessing the thicknesses of the different phases and NMR provides additional information on the co-surfactant distribution given a better contrast between butan-1-ol and water/toluene. With Micro-CT, we found that the final triphasic system is made of 36 % v/v of a bicontinuous phase, 34 % v/v of excess oil phase and 30 % v/v of excess aqueous phase.

Middle bicontinuous phase

As depicted in Figure 4 h., when equilibrium is reached, an homogenous composition of brine, butan-1-ol, SDS and toluene is found in the middle phase. NMR calculations in the middle phase lead to percentages of 47 % v/v of [brine +SDS] (48 % v/v with Micro-CT); 14 % v/v of butan-1-ol ; 39 % v/v of Toluene. The results are coherent with the data provided by Bellocq et al. [START_REF] Bellocq | Three-Dimensional Phase Diagram of the Brine-Toluene-Butanol-Sodium dodecylsulfate System[END_REF] on the same system for a co-surfactant/surfactant ratio of 2 at optimal salinity where 52.5 % w/w of brine and 33.5 % w/w of oil were measured in the middle phase. Similarly, for the composition, no structural gradient was found neither in SAXS nor in DLS (Figure 2). The relaxation time of the bicontinuous middle phase is about 60 µs leading to a diffusion coefficient Dapp = 2.3x10 -11 m²/s. In addition, the SAXS results confirmed the presence of a bicontinuous « sponge-like » structure with a characteristic peak in the intermediate q-range. A correlation length ξ of 13 nm was found and domains periodicity d is 31 nm leading to a deduced amphiphilicity factor fa of -0.75 which indicates that the microemulsion forms a distinct phase [START_REF] Schubert | Small angle neutron scattering near Lifshitz lines: Transition from weakly structured mixtures to microemulsions[END_REF]. DLS and SAXS results are in agreement with the literature for similar systems at optimal salinity [START_REF] Fukumoto | Investigation on Physical Properties and Morphologies of Microemulsions formed with Sodium dodecylbenzenesulfonate, Isobutanol, Brine, and Decane, Using Several Experimental Techniques[END_REF][START_REF] Pouchelon | Light scattering from oil-water interfaces[END_REF][START_REF] Auvray | Concentrated Winsor microemulsions[END_REF][START_REF] Hayes | Observation of a structural gradient in Winsor-III microemulsion systems[END_REF][START_REF] Engelskirchen | Triacylglycerol microemulsions stabilized by alkyl ethoxylate surfactants-A basic study Phase behavior, interfacial tension and microstructure[END_REF].

Composition in the excess oil phase

At equilibrium, no butan-1-ol is detected in the excess oil phase in the limit of detection of our method (Figure 4 h.). Also, the phase does not scatter neither in light nor in the X-ray domain (position F in Figure 7). This confirms that no dispersed object was present in the excess oil phase which is thus composed of toluene.

Composition in the excess aqueous phase

In the excess water phase, butan-1-ol is found at around 3 % v/v with NMR. The value is coherent with Micro-CT where 2 % v/v of other species were detected in the excess aqueous phase. However, this result is not significant regarding the respective accuracy of each method. Nevertheless, it is likely that small quantities of butanol (< 2 % w/w) could remain in the excess aqueous phase as it was previously reported on the same system [START_REF] Bellocq | Three-Dimensional Phase Diagram of the Brine-Toluene-Butanol-Sodium dodecylsulfate System[END_REF][START_REF] Pouchelon | Origin of low interfacial tensions in systems involving microemulsion phases[END_REF]. As the excess oil phase, the aqueous phase does not scatter visible light.

Kinetic study of the microemulsion formation

We illustrated in Figure 3 the different configurations of the sample that were observed in this work during the formation of the microemulsion. The process can be divided into different sections that will be commented below: a. The initial aqueous phase b. The oil/aqueous phase contact c. The original formation of M2 phase d. The evolution of the system into three intermediate phases called "Z3", "M2" and "M1" e. The kinetic aspects leading to the disappearance of the liquid crystal f. The growth of the bicontinuous phase with the fading of M2 and Z3. The mechanisms described here were the same independently of the total volume of the system (~50 µL for SAXS, 3 mL for NMR and 10 mL for DLS and Micro-CT). However, the volume of the system has an impact on its equilibration time. At an oil addition speed of 126.8 mL/hour at ambient temperature, a 3 mL-microemulsion needs 7 days to stabilize while a 10 mL-microemulsion reaches equilibrium in 100 days. Kinetics results for NMR, Micro-CT and DLS are given respectively in Figure 4, Figure 5 and Figure 6. The SAXS acquisitions were performed solely at step 4 (Figure 3) when the three intermediate phases coexisted. The corresponding scattering pictures, acquired in function of height in capillaries, are displayed in Figure 7. SAXS scattering curves for each position in the capillary are shown in Figure 8.

a. Initial aqueous phase

Initially, we added 4 % w/w of SDS, and 10 % v/v of butan-1-ol in 54 g/L of NaCl brine as the aqueous phase (Figure 4 a.). Using DLS and a cumulant analysis, measurements at several values of q confirm a linear relationship between 1/𝜏 and q² in the aqueous phase, characteristic of a diffusive process. The average value of the observed object hydrodynamic diameter is 𝐷 = 40 nm. An SBL-type analysis of the autocorrelation function reveals the presence of polydisperse objects with a size distribution between 𝐷 = 20 nm and 70 nm. However, typical sizes measured for SDS spherical micelles have an hydrodynamic radius 𝐷 of 4-6 nanometers [START_REF] Corti | Quasi-elastic light scattering study of intermicellar interactions in aqueous sodium dodecyl sulfate solutions[END_REF]. As a result, we deduce that the initial aqueous phase at optimal salinity is not mainly made of spherical micelles. It is found in the literature that the addition of salt and alcohol to micellar solution of SDS leads to the formation of rod-like or worms micelles and even surfactant bilayers either vesicles or lamellar phases at higher concentration [START_REF] Oberdisse | Formation spontanée de vésicules dans un système amphiphile chargé: Spontaneous formation of vesicles in a charged amphiphilic system[END_REF][START_REF] Jose | Structural changes of a sodium dodecyl sulfate (SDS) micelle induced by alcohol molecules[END_REF]. Small Angle Neutron Scattering experiments (results to be published) were performed in such phase for wave vectors q between 0.004 and 0.4 Å -1 . It confirmed the presence of cylindrical-shape micelles of 25 nm-length and 2 nm-diameter but did not exclude the presence of larger objects as our lower q-range is limited. The SANS and DLS results are compatible in the sense that the DLS scattering signal accounts for the pseudo-spheres associated with the cylindrical objects.

b. Oil/aqueous phase contact

After phases contact, we noticed that the butan-1-ol progressively migrated at the interface between the aqueous phase and toluene as seen in the NMR compositions profiles (Figure 4 b. and c.). This is characterized by an increase in the brine volume fraction in the aqueous phases between t = 0 to 1 day with our Micro-CT method (see Supplementary Material Part IV). As a confirmation of this effect, the bulk aqueous phase no longer scatters light one day after toluene addition. The upper oily phase also does not scatter light and still corresponds to toluene (Figure 4 b. and c.). We deduced that the surfaceactive species, SDS and butan-1-ol, diffused to the new interface to lower surface tensions when the two initial phases come into contact at t = 0. Visually, the system is still composed of two phases: oil and aqueous phase with an undefined and turbid interface.

c. Original formation of the M2 phase

The third step (Figure 3) is the formation of a first intermediary phase called M2 between the aqueous and oil phases. The situation is illustrated by Figure 4 d. with NMR, Figure 5 with Micro-CT and Figure 6 with DLS. The interfaces of M2 with the oil and aqueous phases are diffuse and opalescent. Macroscopically, the M2 phase is clear and birefringent by observation under crossed polarizers. Surprisingly, at the location of M2, we observed a sharp loss of NMR signal that stems from a drop in T1-relaxation times [34,[START_REF] Gradišek | Study of Liquid Crystals Showing Two Isotropic Phases by 1H NMR Diffusometry and 1H NMR Relaxometry[END_REF], which is displayed as the hatched area (Figure 4 d.). This is due to a drastic reduction of molecules mobility. The local structure is considered as a new species that has a specific relaxation time significantly different from the bulk case. Because this region could not be separately measured, we could not have access to the characteristic T1-relaxation time. Consequently, the complete composition could not be extracted in this zone. We were limited to the fact that M2 is made of a [brine + SDS]-rich phase up to 90 % v/v (Figure 5). Despite its birefringence, from DLS analysis, the M2 phase has a characteristic relaxation time of 71 µs which is the same at the top and bottom of M2 (Figure 6). Due to the evidences of birefringence, its NMR behavior and composition, we deduce that M2 is probably a SDS liquid crystal. Indeed, in the literature, the phase diagram of SDS/water confirms the presence of a semi-crystalline phase at a high concentration of surfactants (~ 30 % w/w) [START_REF] Ekwall | Composition, Properties and Structures of Liquid Crystalline Phases in Systems of Amphiphilic Compounds[END_REF][START_REF] Muzzalupo | Phase Diagram and Dynamic Properties of the Ternary System Water -Sodium dodecylsulfate -Aerosol OT[END_REF][START_REF] Kékicheff | Between cylinders and bilayers[END_REF]. d. Evolution of the system into three intermediate phases: "Z3", "M2" and "M1"

The fourth step is characterized by a more complex intermediate zone (Figure 3). The situation is related to Figure 4 e. in NMR, and between 2 to 8 days in Micro-CT (Figure 5) and DLS (Figure 6). At this stage, the liquid crystal, that was located between the excess oil and water phases, tends to form supernatant and underlying layers with well-defined interfaces. The sample is now composed of excess oil and aqueous phases, and three intermediate phases that are respectively called M1, M2 and Z3 from top to bottom.

We first analyzed the composition of such a system. The excess aqueous and oil phases are mainly composed of brine and toluene respectively. The bicontinuous phase M1 is made of toluene, SDS, brine and a high volume fraction of butan-1-ol (25 % v/v). Butan-1-ol is homogenously and entirely distributed in M1 but small gradients of oil and brine are present (Figure 5 and Figure 4 e.). Regarding M2 phase, significant gradients of composition emerged as seen on the Micro-CT colormap in Figure 5 after about 3 days. The brine volume fraction varies between 60 % v/v and 80 % v/v in M2 from top to bottom of the tube. Finally, Z3 is an homogenous phase containing 96 % v/v of brine given by Micro-CT (see detailed raw profiles in Supplementary Material Part IV). When analyzing the structure and dynamics of the systems, we focused on SAXS and DLS measurements at the formation of the three intermediate phases. The corresponding figures are respectively Figure 7, Figure 8 and Figure 6 at day 7. Excess oil and aqueous phases do not scatter visible light or X-Ray (position F in Figure 7) and are coherent with the found compositions. Concerning M1, we found that it has the same dynamics as the equilibrated microemulsion (Figure 6). We also confirmed that M1 is a bicontinuous microemulsion of which SAXS scattering patterns can be fitted by the Teubner-Strey model (Figure 8). A structural gradient is present in the bicontinuous phase. Indeed, the periodicity of the domains and the correlation length ξ are gradually different from the bottom to the top of the capillary passing from (d=19 nm; ξ=9.5 nm) to (d=22 nm; ξ=11 nm) respectively. We assume that these structural gradients are related to the compositional ones described above. In that sense, a higher concentration of oil would imply larger domains of the bicontinuous network. Additionally, the amphiphilicity factor fa = -0.81 is relatively constant in M1 and smaller than the one measured at equilibrium (fa = -0.75). The structure of the bicontinuous phase, at this time, appears to be more ordered compared to the final one. The decrease of fa in time arises from the increase in domains sizes linked to the incorporation of both water and oil in the bicontinuous phase. In M2, the upper and lower parts of M2 now show distinct values of relaxation times meaning that they are significantly different in terms of dynamics and probably in structure (Figure 6). Under crossedpolarizers, the whole M2 phase is still birefringent and we could clearly identify a change in morphology along the capillary (Figure 7). Indeed, a variety of structures is found in M2 (Figure 7, Position A to C). At the center of M2 (Figure 7 at position A), the azimuthal profile 𝐼(𝜓)corresponds to a vertically oriented liquid crystal phase with a characteristic lamella size (period) of 11 nm (Table 2). The objects have an angle of orientation of 10° and an order parameter of 0.44 (Figure 8), which is characteristic of highly ordered objects. When approaching M1, the objects preserve their sizes but progressively change their orientation (position B, C, D) by rotation until reaching a perpendicular position with an angle of orientation α of 96° (Table 2). The origin of such an orientation remains questionable, it could potentially stem from steric molecular interactions or could be influenced by the initial addition of toluene onto the aqueous phase. The parameter S that characterizes the ordering of the system, i.e. here linked to its anisotropy, decreases when passing from M2 to M1 as described in Table 2. The change in structure goes through a point (B on Figure 7) where both the bicontinuous network and the liquid crystal phase coexist in a similar proportion. At this point, the parameter S is not consistent due to a significant noise in the scattering signals. Above the interface with M2, the structure loses orientation in M1 (position D, E on Figure 7). These are direct evidences of the progressive restructuring of SDS molecules from liquid crystal (in M2) into monolayers (in M1). The sizes of the lamella in M2 are similar to the correlation length ξ of the objects in the bicontinuous phase (~ 10 nm), which is in some cases approximated as the length where the interfacial film is flat [START_REF] Guest | Persistence length measurements in middle phase microemulsions[END_REF][START_REF] Pieruschka | Statistical mechanics of random bicontinous phases[END_REF]. Consequently, it is possible that the bicontinuous network arises from the merging of the sub-units of M2. Regarding the structure type [START_REF] Holmes | Smectic-nematic transition in a lyotropic liquid crystal[END_REF][START_REF] Hendrikx | Anisotropic aggregates of amphiphilic molecules in lyotropic nematic phases[END_REF][START_REF] Hendrikx | Structural relations between lyotropic phases in the vicinity of the nematic phases[END_REF], M2 at position A and B stems from nematic phases as only one and relatively broad scattering peak is detected on the SAXS scattering plots in M2 (Figure 8 at position A, B) [START_REF] Davidson | The measurement of the Nematic Order Parameter by X-ray Scattering Reconsidered[END_REF]. At position C, the peak is sharp and no second order peak is observed. Such a structure is questionable and could be either an electrostatic pseudo-lamellar phase or a smectic phase or an electrostatically swollen sponge phase [START_REF] Davidson | Isotropic, nematic, and lamellar phases in colloidal suspensions of nanosheets[END_REF][START_REF] Jho | Lamellar phase coexistence induced by electrostatic interactions[END_REF] in equilibrium with the sponge phase of position D. Concerning the Z3 phase, it is an opalescent and isotropic zone of relatively constant thickness located between the excess aqueous phase and M2. Z3 is hardly detectable with NMR because it is too close to M2 phase where the signals are strongly impacted by the semi-crystalline structure. In DLS, the simple cumulant inversion method leads to relatively high relaxation times in Z3 phase that are characteristic of large objects in the range of ~100 nm as seen in Figure 6. A polymodal, thus more adapted, Pade inversion analysis was carried out at different angles 𝜃 of 30; 45; 90 and 135°. It follows that the process of relaxation is diffusive and characterized by two diffusion coefficients D ~ 5.10 -11 m².s -1 (𝐷 = 20 𝑛𝑚) and D ~ 5.10 -12 m².s -1 (𝐷 = 200 𝑛𝑚). Because there is one decade between those two coefficients, the result is compatible either with a polydisperse distribution of objects or a bimodal distribution. Regarding the structure, the interpretation of the SAXS data in Z3 is ambiguous as we do not know its dilution regime, and because our lower q-range is limited at 0.01 Å -1 . However, we noted that it could be fitted using a two-population "core-shell structure with thin shell" form factor as detailed in Supplementary Material Part VII. The diameters found are 20 nm (length) and 1.4 nm (radius). DLS and SAXS agree to state that Z3 contains polydisperse core-shell structures but the observation range varied between the two techniques. Further SAXS or SANS experiments would be appreciated to characterize more precisely Z3. We assume that such a phase may appear to ensure the continuity of the thermodynamical properties between the liquid crystal and the excess brine phase. The given analysis leads us to wonder why the system is structured in layers as described above. In the literature, the nucleation of different phases in the vicinity of the original interface has already been observed in various systems and explained by the differences in diffusion coefficients of the species [START_REF] Laughlin | The aqueous phase behavior of surfactants[END_REF] [START_REF] Ricoul | Phase Equilibria and Equation of State of a Mixed Cationic Surfactant-Glycolipid Lamellar System[END_REF]. In our case, the layering is caused by a competition between the migration of the surfactant and the co-surfactant towards toluene. Butan-1-ol is likely to diffuse quickly and totally to the toluene interface immediately after phases contact because of its total solubility and small molecular weight (Mw = 74 g/mol). At the opposite, as a bigger molecule (Mw = 288 g/mol), SDS diffuses with a lower rate to the new interface (see Supplementary Material Part VIII for the reported value of self-diffusion coefficients). The difference between SDS and butan-1-ol diffusivities implies the presence of high composition gradients in the system. Thus, several assemblies of SDS are encountered in function of the height in the system as if the tridimensional phase diagram of the quaternary system was partially covered [START_REF] Bellocq | Three-Dimensional Phase Diagram of the Brine-Toluene-Butanol-Sodium dodecylsulfate System[END_REF][START_REF] Kékicheff | Phase Diagram of Sodium dodecylsulfate-Water System[END_REF][START_REF] Rong | Structure Properties for SDS/ n-C4H9OH/H20 Microemulsion[END_REF]. Indeed, it is known that the molecular exchanges with surfactants occur in the range of the nanosecond to the microsecond [START_REF] Aniansson | Kinetics of step-wise micelle association[END_REF][START_REF] Aniansson | Theory of the kinetics of micellar equilibria and quantitative interpretation of chemical relaxation studies of micellar solutions of ionic surfactants[END_REF]. Thus, the structures observed in M1, M2 and Z3 can locally be considered at equilibrium and correspond to the microstructures of the phase diagram for the compositions that were determined by NMR and Micro-CT. In Figure 9, we indicate in the Toluene/butan-1-ol/6% NaCl Brine phase diagram available in literature [START_REF] Bellocq | Three-Dimensional Phase Diagram of the Brine-Toluene-Butanol-Sodium dodecylsulfate System[END_REF] the points that correspond to t = 0 and the equilibrium as well as the probable location of M2, M1 and Z3 phases. The limit of quantification of our methods prevents us to point out a precise point of the phase diagram. Moreover, we believe that the SDS and butanol composition gradients from bottom to the top of the tube imply that the diffusion path crosses different planes of the butan-1-ol/toluene/Brine phase diagram, each plane corresponds to a specific butan-1-ol:SDS ratio. Because butan-1-ol diffuses towards the toluene supernatant phase, the butan-1-ol:SDS ratio increases from M2 to M1 phase. This conclusion is coherent with literature where "gel" regions were observed for butan-1-ol:SDS ratio lower than 2 (equals to the ratio in the final bicontinuous phase).

Figure 9 The ternary phase diagram n-butanol/toluene/6% NaCl brine for different planes corresponding to a specific butan-1-ol/SDS ratio [START_REF] Bellocq | Three-Dimensional Phase Diagram of the Brine-Toluene-Butanol-Sodium dodecylsulfate System[END_REF]. In our case, n-butanol/SDS ratio is 2 at t=0 and equilibrium. The symbols I, II, III stem from Winsor I, Winsor II, Winsor III region. The dashed zone represents an identified "gel" region in literature and the red zones/arrows represent the probable locations of the M2, M1 and Z3 phases.

e. Kinetic aspects leading to the disappearance of the liquid crystal

In this part, we focused on the evolution over time of the three intermediary phases. This is related to the measurements between day 2 and day 33 in Micro-CT (Figure 5) and DLS (Figure 6). Namely, we extracted the interfaces positions M2/M1 and M1/Oil from Micro-CT profiles and followed their behavior over time. Results are plotted in Figure 10.

Regarding the dynamics evolution of M2 over time, the upper part of M2 shows a maximum in relaxation times and then vanishes when its dynamics becomes comparable to M1 phase (Figure 6). After a few days, the lower part of M2 undergoes the same behavior. We correlate the presence of the maxima in relaxation times to the position of M1/M2 interface over time. The displacement of this maximum indicates the progressive breakage of the liquid crystal phase. Surprisingly, from Figure 10 we found that the displacement of M1/M2 interface is linear in function of time. The flux of matter exchanges through the interface from M2 to M1 is constant at a speed of 4.3×10 -9 m.s -1 . Curves plots are displayed in Supplementary Material Part V for more details.

Concerning M1, its thickness increases by progressive solubilization of oil and brine by conserving a mean water/oil ratio around 50/50 v/v (Figure 5). Also, the displacement of the oil/M1 interface is proportional to √𝑡 and thus follows Fick's law with a characteristic diffusion coefficient of 1.2 x10 -11 m².s -1 as seen in Figure 10. This is in agreement with Broens et al. study in flowing conditions [START_REF] Broens | Emulsification kinetics during quasi-miscible flow in dead-end pores[END_REF]. It means that the penetration of oil is driven by the composition gradients present in the bicontinuous phase. In addition, the DLS monitoring highlights that M1 reaches a minimum in relaxation times at the total disappearance of M2 (Figure 6). An increase of the ordering of the bicontinuous network (i.e decrease of fa or increase of ξ/d) could explain the decrease in relaxation times. Indeed, the correlation between ξ/d and the relaxation time τ has already been found in previous work by Fukumoto et al. [START_REF] Unsal | Imaging of compositional gradients during in situ emulsification using X-ray micro-tomography[END_REF].

In that sense, this result would particularly support the static SAXS analysis where the phase M1 was found to be more ordered than the final bicontinuous phase.

Overall, given the evidences of the dynamic and static analysis, we confirmed that the bicontinuous phase M1 is growing towards the excess aqueous interface by progressive merging with the liquid crystal M2. The liquid crystal M2 then acts like a reservoir of surfactants for the building of the bicontinuous phase M1 (i.e. the future middle phase of the equilibrated system). f. The growth of the final bicontinuous phase

In the last step of the process (Figure 3), the growth of M1 finally leads to the total consumption of the liquid crystal M2. The disappearance of the liquid crystal M2 causes the inherent fading of the Z3 phase (Figure 5). This provides another source of surfactants to the bicontinuous phase M1. As a consequence, it enables the additional solubilization of toluene and water that leads to a thickening of the middle phase until equilibrium is reached (Figure 4. f., g., h. and from day 29 in Figure 5). We noticed a symmetrical growth of the middle phase along the tube leading to a centered middle phase in the final system as expected for a microemulsion at optimal salinity [START_REF] Huh | Interfacial tensions and solubilizing ability of a microemulsion phase that coexists with oil and brine[END_REF]. Structural and compositional properties of the final system were previously commented in 3.1.

Conclusions

Applied to Enhanced Oil Recovery, previous kinetic studies on microemulsion formation were performed in dynamic conditions to better understand the phases behaviors of such systems in flowing conditions [START_REF] Unsal | Pore Scale Dynamics of Microemulsion Formation[END_REF][START_REF] Broens | Emulsification kinetics during quasi-miscible flow in dead-end pores[END_REF][START_REF] Unsal | Imaging of compositional gradients during in situ emulsification using X-ray micro-tomography[END_REF][START_REF] Tagavifar | Spontaneous and Flow-Driven Interfacial Phase Change[END_REF][START_REF] Shahidzadeh | Spontaneous emulsification relation to microemulsion phase behaviour[END_REF][START_REF] Alzahid | Flow regimes during surfactant flooding[END_REF][START_REF] Yu | Formation and Flow Behaviors of in Situ Emulsions in Heavy Oil Reservoirs[END_REF], especially at the vicinity of dead-end pores. However, the information extracted by such an investigation is limited due to complexity of model rocks or non-applicability of quantitative advanced techniques. To extend the knowledge on the in-situ mechanisms, we monitored in bulk the formation of a bicontinuous microemulsion from oil/water contact to equilibrium. We succeeded in using DLS, Micro-CT, NMR [START_REF] Herrera | Quantification of microemulsion systems using low-field T1-weighted imaging[END_REF] and SAXS techniques in a non-destructive way to quantitatively analyze the properties of the system over time. First, the techniques were implemented on an equilibrated sample to assess the methods and to have access to the structure and compositions of the final system. A good agreement was found with the literature on a similar system [START_REF] Bellocq | Three-Dimensional Phase Diagram of the Brine-Toluene-Butanol-Sodium dodecylsulfate System[END_REF][START_REF] Pouchelon | Light scattering from oil-water interfaces[END_REF][START_REF] Auvray | Concentrated Winsor microemulsions[END_REF][START_REF] Pouchelon | Origin of low interfacial tensions in systems involving microemulsion phases[END_REF] which made it possible to carry out the kinetic study. The work shows that the formation of the model microemulsion involves the presence of transitional layers in which structures/dynamics properties and compositions are correlated according to the position in the phase diagram [START_REF] Bellocq | Three-Dimensional Phase Diagram of the Brine-Toluene-Butanol-Sodium dodecylsulfate System[END_REF]. We were able to identify different steps of the microemulsion formation. The initial aqueous phase contains cylindrical objects that migrate at the interface when in contact with oil. Then, an intermediate liquid crystal appears and causes the formation of the supernatant and underlying layers, respectively a bicontinuous phase and a dispersed phase. The final bicontinuous phase of the Winsor III system stems from the dismantling of the liquid crystal phase. The latter acts as a reservoir of surfactants by providing a constant matter flow toward the bicontinuous phase until total consumption. The sub-units of the liquid crystal are about 10 nm and can be compared to bricks that build up the bicontinuous network. Several types of liquid crystals that have different orientation angles were observed and should be studied in a more detailed SAXS analysis. Toward the excess oil phase, the penetration front [bicontinuous phase/excess oil] follows the Fick's law with a diffusion coefficient that is comparable to those observed under flow conditions in dead-end pores [START_REF] Broens | Emulsification kinetics during quasi-miscible flow in dead-end pores[END_REF]. Toward the water excess phase, the liquid crystal coexits with a dispersed phase that potentially contains large vesicles. Rheological studies would be appreciated to assess the viscosity and interfacial properties of the found liquid crystal phase -provided that the liquid crystal structure is preserved when sampled. Applied to Enhanced Oil Recovery, this would give further information on the flowing behavior of such a phase in-situ in the reservoir. Indeed, the formation of a liquid crystal could potentially lead to a reservoir plugging and a decrease in oil solubilization performances at larger scale. Also, this study should provide interesting routes to understand further the mechanisms of formation of a Winsor I or II system where spherical structures are found at equilibrium. One can wonder if they also arise from a crystalline phase. It could eventually also be extended to non-ionic surfactants, i.e. when no co-surfactant is added to the formulation, to stress their roles in the mechanism of formation and have a better insight on the nature of phase transitions occurring in-situ.

Figure 1

 1 Figure 1 Photograph of the SDS/butan-1-ol/NaCl brine/toluene microemulsion at optimal formulation at equilibrium and an intermediate state

Figure 2

 2 Figure 2 Characteristics of the equilibrated middle phase from SAXS and DLS measurements. Left: SAXS in absolute intensity (cm -1 ) versus the scattering vector q (nm -1 ) at the bottom and top in the middle phase fitted with Teubner-Strey model ; right: DLS Normalized Autocorrelation function G (a.u) in function of time (µs) at the bottom and top in the middle phase.

Figure 3

 3 Figure 3 Schematic description of the phases behavior in the system during the formation of the microemulsion -from the Toluene addition to the equilibrium

Figure 4

 4 Figure 4 Composition profiles after Toluene addition onto the aqueous phase until equilibrium is reached -Total volume 3 mL-calculated from NMR 1D-relaxation profiles. The dashed zones correspond to a region where quantification is not possible due to the presence of the liquid crystal phase.

Figure 5

 5 Figure 5 Brine volume fraction mapping in function of height in the tube and time (Total volume = 10 mL) -calculated by Micro-CT from the corresponding radiographies of the tube in time in level of grey, displayed above. Vertical lines show the times where the acquisitions were carried out. White lines represent the composition contours as guidelines. The dashed horizontal white line is indicative of the interface of Z3 with the excess aqueous phase as it cannot be highlighted with the chosen color map. (See in Supplementary Material Part IV that details a better representation of Z3).

Figure 6

 6 Figure 6 Evolution of relaxation times during the formation of the microemulsion in the intermediate phases from top (M1) to bottom (Z3) (Total Volume = 10 mL) -Extracted from a classical cumulant method from Dynamic Light Scattering measurements at 𝜃= 90°-Vertical lines represent the acquired measurements times.

Figure 7

 7 Figure 7 Schematic description of the SAXS experiment on a non-equilibrated microemulsion. (left) Location of the measured scattering points and corresponding transmissions in function of height z in the tube -(middle) Observation of the M2 phase in the SAXS capillary under crossed-polarized light. -(right) The scattering acquisitions that are represented at each measurement point in the sample.

Figure 8

 8 Figure 8 Scattering profiles according to the position in the capillary introduced in Figure 7. Scattering curves in M1 (D, E) were fitted by the Teubner-Strey model; Scattering curves in M2 (A, B, C, (D)) with Maier-Saupe fits. Scattering pattern for point G (divided by 5 for better representation) was fitted with a two-populations core-shell form factor (see in Supplementary Material Part VII).

Figure 10

 10 Figure 10 Relative position X of M2/M1 and M1/Oil interfaces (normalized to initial position) in function of time and fits associated.

Table 1

 1

	% volume	Phase	Brine + SDS	Others
		Excess oil	(2)		98
	Micro-CT ( +/-2 % v/v)	Bicontinuous	48		52
		Excess aqueous	98		(2)
			Brine	SDS	butan-1-ol	Toluene
		Excess oil				100
	NMR ( +/-4 % v/v)	Bicontinuous	42	5*	14	39
		Excess aqueous	97		(3)
	() Not significant
	* Determined by hypothesis

Composition of the equilibrated SDS/butan-1-ol/NaCl brine/Toluene microemulsion at optimal formulation -from NMR and Micro-CT

Table 2

 2 Parameters extracted from the Maier-Saupe fits at position A, B, C, D in M2/M1 region introduced in Figure 7. Raw data are available in Supplementary material Part VI

	Position	Order parameter S	Orientation angle α (°)	Periodicity (nm)
	A	0.44	10	11
	B	(0.04)	51	11
	C	0.14	77	11
	D	0.063	96	11