Tarik Chakkour

Emmanuel Frénod

Inverse problem and concentration method of a continuous-in-time financial model *

In a continuous-in-time model there is the important financial quantity called Loan which can not be determined directly in terms of Algebraic Spending but has a major impact on the financial strategy. In this paper, we use a mathematical framework to discuss an inverse problem of determining the implied Loan Measure from Algebraic Spending Measure when it is possible. In addition, we build a numerical method to concentrate a measure as a sum of Dirac masses.

Introduction

The continuous-in-time financial model, first introduced by Sundaresan in [START_REF] Suresh | Continuous-time methods in finance: A review and an assessment[END_REF], constitutes a powerful tool for studying the development of continuous-in-time methods in finance. These include computational and estimation methods to test and implement continuous-time models. We refer to papers [START_REF] Chiarella | The calibration of stock option pricing models using inverse problem methodology[END_REF][START_REF] Robert | Continuous-Time Finance[END_REF], which are dealing with continuous-in-time financial model.

At the same time, inverse problems began to be studied and applied systematically in various branches of modern science as in finance. Studying problems can be classified as inverse, and they are among the most complicated ones (since they are unstable and usually nonlinear). The reason is that solutions to inverse problems describe important properties of variable under study such as density or measure and the integrability property. Many studies have used an inverse problem in finance. For example the theory of inverse problems via Tikhonov regularisation in [START_REF] Egger | Tikhonov regularization applied to the inverse problem of option pricing: convergence analysis and rates[END_REF] is used in order to analyse the procedure for determining volatility in the model of Black and Scholes. H. Egger, H.W. Engl [START_REF] Frénod | A continuous-in-time financial model[END_REF] study the inverse problem of option pricing under consideration and they prove a stable solution of the inverse problem and its convergence.

In recent paper [START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF] we have also constructed a continuous-in-time model which is designed to be used for the finances of public institutions and not for the financial market. This model permit to set out annual and multiyear budgets for any organization and describes working of loan and repayment in order to forecast its future financial plans. We build this model because of the default of a discrete model which is using tables. Indeed, we justified in the introduction of [START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF] why that we created this model. This model relies on measure theory and uses the mathematical tools such convolution and integration, etc. In addition, we check the mathematical consistency of this model in [START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF]. Frénod & Safa [START_REF] Frénod | Two optimization problems using a continuous-in-time financial model[END_REF][START_REF] Robert C Merton | Theory of finance from the perspective of continuous time[END_REF][START_REF] Frńod | Continuous-in-time financial model for public communities[END_REF] use this model in the framework of the control theory in order to set out an optimal control method for the strategy elaboration phase to better adjust the project implementation.

In a continuous-in-time model an Algebraic Spending takes into account the spending and the income to satisfy the current needs. There exists a mathematical operator that links between Algebraic Spending and Loan. If this operator is not invertible, solutions of the posed inversion problem can be brought.

The aim of the present paper is to study the inverse problem with giving a rigorous mathematical formulation to determine Loan from the Algebraic Spending and to build a concentration method in order to concentrate a measure at positions. The paper is organized as follows. We begin by describing in section 2 a continuous-in-time model. In section 3, we study the inverse problem that link between financial quantities in the model. We analyse it on measure space in order to show its consistency. In section 4, we introduce a numerical approach which is allowing to concentrate any measure as a combinaison of Dirac measures.

Continuous-in-time financial model

In this section we recall the variables used in the model and introduced in [START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF] in order to describe links existing between them and their meanings. These variables are defined on the set of Radon measures.

The Radon measure space

We build the continuous-in-time financial model on the Radon measure space M([t I , Θ max ]) which is a continuous and linear form acting on continuous functions space

C o c ([t I , Θ max ]) defined over a time interval [t I , Θ max ]. The usual norm on M([t I , Θ max ]) is: µ M((tI,Θmax)) = sup ψ∈C o c ([tI,Θmax]),ψ =0 µ, ψ ψ L ∞ ([tI,Θmax]) , (1)
where

L ∞ ([tI,Θmax]) is the usual norm on C o c ([t I , Θ max ]) defined as: ψ L ∞ ([tI,Θmax]) = sup t∈[tI,Θmax] |ψ(t)| . ( 2
)
We notice that we build the model on a Banach space M([t I , Θ max ]) in order to use the convolution and Fourier Transform.

Financial quantities of the model

We explained in [START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF] the financial quantities that are involved in the models and the relations between them. For instance the Loan Measure κE and the Repayment Measure ρK are connected by a convolution operator:

ρK = κE γ, (3)
where the Repayment Pattern Measure γ is a non-negative measure with total mass which equals 1, γ has to satisfies:

+∞ -∞ γ = 1. ( 4
)
We denote by K RD the Debt Field, given as the solution to the following Ordinary Differential Equation:

dK RD = κE -ρK -ρI K , (5)
where measure ρI K is the repayment of the Current Debt K RD at initial time t I . It is called initial Debt Repayment Scheme. Initial condition for equation ( 5) is given by:

+∞ tI ρI K = K RD (t I ). ( 6
)
The solution of this ODE is expressed:

K RD (t) = K RD (t I ) + t tI κE - t tI ρK - t tI ρI K = t tI κE - t tI ρK + +∞ t ρI K . ( 7
)
We denote by ρI the Interest Payment Measure defined by:

ρI = αK RD (t)dt, ( 8
)
where α is the floating rate interest. Algebraic Spending Measure σ is defined such that the difference between spendings and incomes required to satisfy the current needs. Algebraic Spending Measure σ takes into account the spending and the income and it is given by:

σ = σg + β, ( 9
)
where σg is the Current Spending Measure and β is the Measure of Isolated Spending. The balance equation between incomes and spendings is expressed by the following equality:

κE = σ + ρK + ρI K + ρI . ( 10
)

Inverse problem of the model

In this section we inverse operators given below which are useful of the model. We study this inversibility in measure space. For that, we consider a positive number Θ γ such that Θ γ < Θ max -t I and we set the Repayment Pattern Measure γ such that:

γ ∈ M([0, Θ γ ]). ( 11
)
In order that Repayment Measure ρK given by a convolution operator in relation [START_REF] Frénod | A continuous-in-time financial model[END_REF]

σ = κE -κE γ -α κE -κE γ, 1 |[tI,t[ dt -α ρI K , 1 |[t,Θmax[ dt -ρI K . ( 13
)
Let L be a linear operator defined from M([t I , Θ max -Θ γ ]) to M([t I , Θ max ]) acting on Loan Measure κE by

L[κ E ] = κE -κE γ -α κE -κE γ, 1 |[tI,t[ dt, (14)
where κE -κE γ,

σ = L[κ E ] + D[ρ I K ]. ( 18
)
We define a density measure ẽα with its density e α is exponential function of coefficient α, i.e ẽα = e α (t)dt, where e α (t) = e αt .

Lemma 3.2. If Loan

Measure κE is in M([t I , Θ max -Θ γ ]
) and if γ satisfies relation (11), then we have the following equality:

(1 -F(γ))F(κ E ) = F L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[ , ( 19
)
where F stands for the Fourier Transform Operator, where operator L is defined by (14).

Proof. We have the following equality:

κE -κE γ, 1 |[tI,t[ dt = ẽα κE -κE γ -α κE -κE γ, 1 |[tI,s[ ds, e -α |[tI,t[ . ( 20
)
Indeed, relation (20) is obtained using the integration by parts:

κE -κE γ, αe α|[t I,s[ ds, e -α |[tI,t[ dt = ẽα κE -κE γ, e -α |[tI,t[ -κE -κE γ, 1 |[tI,t[ dt. ( 21
)
Using definition (14) of operator L, equality in (20) is multiplied by α to give:

κE -κE γ, α1 |[tI,t[ dt = αẽ α L[κ E ], e -α |[tI,t[ . ( 22
) Replacing κE -κE γ, α1 |[tI,t[ dt in relation (22) by κE -κE γ -L[κ E ],
we obtain the following equality:

κE -κE γ = L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[ . ( 23
∃ > 0, 1 1 -F(γ) ∈ L ∞ ((-∞, -) ∪ ( , +∞)), (24)
and if Initial Debt Repayment Measure ρI K and Algebraic Spending Measure σ are in M([t I , Θ max ]) and satisfy the following equality:

σ -D[ρ I K ] + αẽ α σ -D[ρ I K ], e -α |[tI,t[ , 1 |[tI,Θmax[ = 0, ( 25
)
then, there exists an unique Loan Measure κE in M([t I , Θ max -Θ γ ]) (if it exists) given by following formula:

κE = F -1 F σ -D[ρ I K ] + αẽ α σ -D[ρ I K ], e -α |[tI,t[ 1 -F(γ) , ( 26
)
such that (18) holds.

In order to prove this theorem, we need the following lemmas:

Lemma 3.4. If
Ker(L) = {κ E ∈ M([t I , Θ max -Θ γ ]), κE -κE γ = 0}. ( 28
)
If Repayment Pattern Measure γ is the Dirac mass δ 0 , then, relation (28) gives that Ker(L) is equal to M([t I , Θ max -Θ γ ]). In fact, δ 0 is the identity element for the convolution operator. Thus linear operator L is not a one-to-one map if γ = δ 0 . We distinguish four cases to prove the lemma:

First case:

In the first case, we will show that if Repayment Pattern Measure γ is density measure satisfying relation (11), then, linear operator L given in relation ( 14) is a one-to-one map.

For any positive integer n, we define V n the vector space of density measures such its densities are polynomial functions in x of degree n and restricted to interval (t I , Θ max -Θ γ ). By density of space V m for any integer m in M([t I , Θ max -Θ γ ]), we will show that for the following Repayment Pattern Measure γ ∈ V m i=m i=0

α i t i i=m i=0 α i Θ i+1 γ i + 1 × 1 [0,Θγ ] (t)dt, ( 29
) that Ker(L) ∩ V n = {0}. (30)
In order to show relation (30), we will show that the coefficients of density measure κE are zero for each κE in Ker(L)

∩ V n . If κE ∈ Ker(L) ∩ V n , then, for any x ∈ (t I , Θ max -Θ γ ), κE = i=n i=0 c i x i × 1 [tI,Θmax-Θγ ] (x)dx, (31)
where (c i ) 0≤i≤n are its coefficients and

κE -κE γ = 0. ( 32
)
As γ is given in relation (29), expression κE γ can be simplified with using relation (31

) κE γ = dx i=m i=0 α i Θ i+1 γ i + 1 × i=n i=0 c i j=m j=0 α j [tI,Θmax-Θγ ]∩[x-Θγ ,x] y i (x -y) j dy. ( 33
) As (x -y) j = p=j p=0 C p j (x) j-p (-y) p , ( 34
) expression (33) gives κE γ = dx i=m i=0 α i Θ i+1 γ i + 1 × i=n i=0 c i j=m j=0 α j p=j p=0 C p j (-1) p (x) j-p × [tI,Θmax-Θγ ]∩[x-Θγ ,x] y i+p dy . (35)
Now, on the one hand, we show that coefficients (c i ) 0≤i≤n are zero, when t I + Θ γ ≤ Θ max -Θ γ . In this case, if x is chosen outside the support of κE such that t

I + Θ γ ≤ Θ max -Θ γ < x < Θ max , then, κE is null. Thus, the intersection of intervals [t I , Θ max -Θ γ ] and [x -Θ γ , x] is interval [Θ max -Θ γ , x -Θ γ ].
Then, relations (32) and ( 35) are used with the fact that Loan Measure κE is null in order to obtain the following equality:

i=n i=0 c i j=m j=0 α j p=j p=0 C p j (-1) p (x) j-p i + p + 1 × (Θ max -Θ γ ) i+p+1 -(x -Θ γ ) i+p+1 = 0. ( 36
)
From relation (36), we get a nul polynom on interval (Θ max -Θ γ , Θ max ) which has an non empty interior. Consequently, we showed that coefficients (c i ) 0≤i≤n are zero, when

t I + Θ γ ≤ Θ max -Θ γ .
On other hand, we show that coefficients (c i ) 0≤i≤n are zero, when Θ max -

Θ γ < t I + Θ γ . Then, if x is chosen outside the support of κE such that Θ max -Θ γ < t I + Θ γ < x < Θ max , then, κE is null. Thus, the intersection of intervals [t I , Θ max -Θ γ ] and [x -Θ γ , x] is interval [Θ max -Θ γ , x -Θ γ ]
. We obtain the same null polynom in relation (36) on interval (t I + Θ γ , Θ max ) which has an non empty interior.

Consequently, coefficients (c i ) 0≤i≤n are zero. From this, we conclude that (30) is true for any integer n. The proof of the first case is achieved.

Second case:

In the second case, we will show that if Repayment Pattern Measure γ is a density measure satisfying relation (11), then, linear operator L given by relation ( 14) is a one-to-one map for any Loan Measure κE which is supposed a concentrated measure (as like as a combination of Dirac measures or a combination of Dirac measures multiplied by a continuons functions). In fact, since the convolution of γ and κE is concentrated measure, then, relation (27) gives that a density measure is equal to a concentrated measure. This equality is possible, if and only if, Loan Measure κE is zero. Consequently, the proof of the second case is achieved.

Third case:

We show in the third case that if Repayment Pattern Measure γ is a concentrated measure satisfying relation (11), i.e: γ = j=m j=1 e j δ aj .

(37)

Thus, for any Loan Measure κE which is supposed a density measure, we have:

κ E (x) = j=m j=1 e j κ E (x + a j ). ( 38
)
Then Loan Density κ E is null. From this, linear operator L given by relation ( 14) is a one-to-one map for any Loan Density Measure κE .

Fourth case:

In the fourth case, we set the Repayment Pattern Measure γ as a combination of m Dirac measures located at points (a j ) 1≤j≤m that is different from δ 0 and satisfying relation (11) such that: γ = j=m j=1 e j δ aj .

(39)

We will show that the coefficients of concentrated measure κE are zero for any Loan Measure κE that can be written in the following form κE = i=n i=1 f i δ bi , where points b i are in [t I , Θ max -Θ γ ] for i from 1

to n. According to (27), we have:

i=n i=0 f i × (δ bi - j=m j=0 e j δ bi+aj ) = 0. ( 40
)
In order to show that coefficients (f i ) 0≤i≤n are zero, we distinguish two situations of points (a j ) 1≤j≤m . The first situation is when points (a j ) 1≤j≤m are not nulls. Thus, since points (e j ) 0≤j≤m are fixed points, relation (40) proves that coefficients (f i ) 0≤i≤n are zero. The second situation is when there exists

a k = 0 for 1 ≤ k ≤ m, then relation (40) gives i=n i=0 f i × (1 -e k )δ bi - j=m j=0,j =k e j δ bi+aj = 0. (41)
Since Repayment Pattern Measure γ is different from δ 0 , relation (41) proves that coefficients (f i ) 0≤i≤n are zero.

Consequently, the study of these four cases is achieving the proof of the lemma.

Lemma 3.5. Assuming that:

γ, t1 |[tI,Θmax[ = 0, ( 42
)
is achieved and if operator L given by relation (14) satisfies:

L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[ , 1 |[tI,Θmax[ = 0, (43)
then, F(κ E ) ∈ L ∞ (R) and is such that

lim ξ→0 F(κ E )(ξ) = - L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[ , t1 |[tI,Θmax[ γ, t1 |[tI,Θmax[ . ( 44
)
If L does not satisfy the equality in relation (43), then, F(κ E ) has an infinite limit in 0.

Proof. As L[κ E ] ∈ M([t I , Θ max ]
), we get:

L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[ ∈ M([t I , Θ max ]). ( 45
) Function F L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[
is expanded in 0 until the order 1 with using the Taylor expansion of exponential function e -iyξ in 0 until the order 1 to get:

F L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[ (ξ) = L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[ , 1 |[tI,Θmax[ -iξ L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[ , t1 |[tI,Θmax[ + O(ξ 2 ) (46)
Since operator L satisfies equality in relation (43), relation (46) gives:

F L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[ (ξ) = -iξ L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[ , t1 |[tI,Θmax[ + O(ξ 2 ) ( 47
)
According to (11), function 1 -F(γ) is Taylor expanded in 0 until the order 1 to obtain:

1 -F(γ)(ξ) = iξ γ, t1 |[tI,Θmax[ + O(ξ 2 ). ( 48
)
According to relations (19), ( 47) and (48), we get equality (44).

As L[κ E ] and γ are integrables on their supports, the numerator and the denominator of right term in relation (44) are finites. From this and according to relation (42), function F(κ E ) admits a finite limit at 0. It is concluded

F(κ E ) is in space L ∞ (R).
On other hand, if equality in relation ( 43) is not satisfied, then, according to relations (19), ( 46) and (48), F(κ E ) has an infinite limit at 0. From this, the proof of lemma is achieved. Lemma 3.6. Under assumptions (11), ( 24), (43) and if operator L given by relation (14) satisfies:

F(L[κ E ]) ∈ L 2 ([t I , Θ max ]), ( 49
)
then, F(κ E ) ∈ L 2 (R).

Proof. According to (49), we get:

F L[κ E ] + αẽ α L[κ E ], e -α |[tI,t[ ∈ L 2 (R). (50)
We notice that for any positive , the product of functions in L ∞ ((-∞, -) ∪ ( , +∞)) and in L 2 (R) is a function in L 2 ((-∞, -) ∪ ( , +∞)). From relations (19), (50) and according to the assumption of γ given in relation (24), we have:

∃ > 0, F(κ E ) ∈ L 2 ((-∞, -) ∪ ( , +∞)). ( 51
)
According to Lemma 3.5, F(κ E ) admits a finite limit at point ξ = 0. Consequently, we have the following result:

∃ > 0, F(κ E ) ∈ L 2 ([-, ]). ( 52
)
From this and using relation (51), the lemma is proved.

Proof of Theorem 3.3. We notice that Repayment Pattern Measure γ equals to Dirac mass δ 0 doesn't satisfy relation (24). From Lemma 3.4, we obtain the uniqueness of κE , achieving the proof of the theorem.

A concentration method of the model

We developped a new approach of concentrating a measure to enrich the model. This approach is allowing to generate the result measure which is called the concentrated measure. Indeed, concentrate Repayment Measure ρK can be interpreted from the financial slant as follow. It means that Repayment amount presented by ρK is associated with localized actions or payments at differents instants in order to express how these payments are distributed at these instants. Thus, the total repayment amount is equal to the sum of these amounts. This equality defines a principle of conservation that we explain it below.

For any integer n strictly greater than 1, we define J the set of integers between 1 and n, i.e J = [ [1; n]]. To concentrate a given measure m, we assume that its integral over R is finite

+∞ -∞ m < ∞. ( 53
)
The concentrated measure of measure m is defined as a sum of Dirac masses (M i ) i∈J located in given positions. The basic idea of concentrating the measure m is to determine the adequate partition of unity which is defined by a set of nonnegative compactly supported functions (φ i ) i∈J . Afterwards, measures defined by a multiplication of measure m by functions (φ i ) i∈J are integrated in order to compute Dirac masses (M i ) i∈J .

To construct this partition of unity, we need three quantities. The first one is the set of intervals ([d i , f i ]) i∈J , where d i < f i for each i ∈ J. In order that the intersection of intervals [d i , f i ] and [d i+1 , f i+1 ] is non empty for each i ∈ J \ {n}, we state the following assumption:

∀i ∈ J \ {n}, d i+1 < f i . ( 54
)
The second one is the fuzzy index F I which influences on calculated Dirac masses

(M i ) i∈J .
i , f i ] for each i ∈ J ∀i ∈ J, Supp(φ i ) = [d i , f i ]. ( 55
)
This allows us to define Dirac masses (M i ) i∈J as follows:

∀i ∈ J, M i = fi di m × φ i , ( 56
)
where, a measure m × φ i is the product of measure m by function φ i . Now, we show the principle of conservation of masses which states that the sum of Dirac masses (M i ) i∈J and the total mass of measure m presented by its integral are equals. Since functions (φ i ) i∈J form a partition of unity, we have i∈J φ i = 1 which is multiplied by measure m to get:

m = i∈J m × φ i . ( 57
)
By integrating over R measure m given by relation (57) and by permuting a finite sum and finite integral, we get the following relation:

+∞ -∞ m = i∈J M i , ( 58
)
which expresses that the total mass +∞ -∞ m is the sum of the Dirac masses (M i ) i∈J . It is concluded that the principle of conservation is proved.

Construction of the partition of unity

It is difficult to build explicitly the partition of unity in our model. Consequently, this partition of unity is determined with approximated numerical method in two steps. The first step is an initialization and the second step is an iteration. We explain these steps as follows:

Initialization step

The aim of the initialization step is to initialize positives functions (φ 0 i ) i∈J defined on R which are not necessary a partition of unity. We notice that the expression of functions (φ 0 i ) i∈J depends on the fuzzy index F I and on its support which is given by following relation:

∀i ∈ J, Supp(φ 0 i ) = [d i , f i ]. (59)
We determine functions (φ 0 i ) i∈J in three cases of the fuzzy index F I which are F I = 0, F I = 1 and F I ≥ 2.

If F I = 0, function φ 0 i (see Figure 1) is a constant piecewise function, which is 1 over

[d i , f i ] for each i ∈ J. ∀i ∈ J, φ 0 i (x) = 1 [di,fi] (x). (60)
x di fi

• • φ 0 i 1 0 y Figure 1: Definition of functions (φ 0 i ) i∈J when F I = 0
If F I = 1, we define the point c i as the mid-point of the line segment [d i , f i ]. Then, function φ 0 i (see Figure 2) is an affine piecewise function on intervals [d i , c i ] and [c i , f i ] for each i ∈ J, which is defined as follows. The restriction of function φ 0 i on interval [d i , c i ] is an affine function which is 0 at point d i and is 1 at point c i for each i ∈ J. The restriction of function

φ 0 i on interval [c i , f i ] is an affine function which is 1 at point c i and is 0 at point f i for each i ∈ J. ∀i ∈ J, φ 0 i (x) = x -d i c i -d i 1 [di,ci[ (x) + x -f i c i -f i 1 [ci,fi] (x). (61)
If F I ≥ 2, function φ 0 i (see Figure 3) is an affine piecewise function on intervals [d i , c i ] and [c i , f i ] for each i ∈ J, where the points d i and f i are to be determined respectively in terms of values (d i ) i∈J and (f i ) i∈J and where the point c i is the mid-point of the line segment [d i , f i ]. Now we will determine the points d i and f i for each i ∈ J as follows. We define each value d i by a value

dj x di ci fi • • • φ 0 i 1 0 y Figure 2: Definition of functions (φ 0 i ) i∈J when F I = 1 ∀i ∈ J, d i = { d j 2 F I , min dj (d i -d j ) > 0} ∪ { d i 2 F I }. (62)
Similarly for the point d i , the point f i is determined by a value fj 2 F I , where point f j is the right neighborhood of the point f i such that a distance f j -f i is strictly positive and is minimized. If there is any point strictly superior than f i , we take f j = f i . Formally, points (f i ) i∈J are defined as follows:

∀i ∈ J, f i = { f j 2 F I , min fj (f j -f i ) > 0} ∪ { f i 2 F I }. (63)
Since points (d i ) i∈J and (f i ) i∈J are given respectively by relations ( 62) and (63), we determine function φ 0 i on interval [d i , f i ] for each i ∈ J when F I ≥ 2 as follows. The restriction of function φ 0 i on interval [d i , c i ] is an affine function which is 0 at point d i and is 1 at point c i for each i ∈ J. The restriction of function φ 0 i on interval [c i , f i ] is an affine function which is 1 at point c i and is 0 at point f i for each i ∈ J. Formally, functions (φ 0 i ) i∈J are defined as follows:

∀i ∈ J, φ 0 i (x) = x -d i c i -d i 1 [d i ,c i [ (x) + x -f i c i -f i 1 [c i ,f i ] (x). ( 64
)
x

d i c i f i • • • φ 0 i 1 0 y Figure 3: Definition of functions (φ 0 i ) i∈J when F I ≥ 2

Iteration step

The aim of the iteration step is to compute functions (φ i ) i∈J which form a partition of unity in terms of initialized functions (φ 0 i ) i∈J given by relation (60) if F I = 0 or by relation (61) if F I = 1 or by relation (64) if F I ≥ 2. Since (CL i ) i∈J is an optional parameter in our modeling, functions (φ i ) i∈J are determined in two cases. The first case is when we don't utilize (CL i ) i∈J . The second case is when we utilize (CL i ) i∈J . In the first case, functions (φ i ) i∈J are defined as follows:

∀i ∈ J, ∀x ∈ i∈J [d i , f i ], φ i (x) = φ i 0 (x) i∈J φ i 0 (x) . ( 65
)
In the second case, functions (φ i ) i∈J are defined as follows:

∀i ∈ J, ∀x ∈ i∈J [d i , f i ], φ i (x) = CL i × φ i 0 (x) i∈J CL i × φ i 0 (x) . ( 66
)
Now, we can verify that functions (φ i ) i∈J given by relation (65) or by relation (66) form a partition of unity. Indeed, functions (φ i ) i∈J are positives because of the positivity of functions (φ 0 i ) i∈J and of values (CL i ) i∈J and the sum of functions (φ i ) i∈J is 1 satisfying following equality:

i∈J φ i (x) = 1.
(67)

Examples of the partition of unity

In order to illustrate the construction of the partition of unity developped in subsection 4.1, we draw examples of functions φ 1 and φ 2 in Figures 4 and5 which are defined on R with support being respectively [0, 6] and [START_REF] Egger | Tikhonov regularization applied to the inverse problem of option pricing: convergence analysis and rates[END_REF][START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF], i.e d 1 = 0, f 1 = 6, d 2 = 2 and f 2 = 4. In addition, since the fuzzy index F I is equal to 0 in simulation of Figure 4, functions φ 1 and φ 2 are piecewise constants functions and are computed with relation (65) or relation (66).

The left picture of Figure 4 shows the partition of unity formed by functions φ 1 and φ 2 which are computed with relation (65). The top diagram presents function φ 1 which is equal to 1 2 between instants 2 and 4 and to 1 between instants 0 and 2 and between instants 4 and 6 and to 0 elsewhere. The middle diagram presents function φ 2 which is equal to 1 2 between instants 2 and 4 and to 0 elsewhere. The bottom diagram shows that function φ 1 + φ 2 is equal to 1 between instants 0 and 6. This example illustrates that formula (65) computes correctly the partition of unity.

The right picture of Figure 4 shows the partition of unity formed by functions φ 1 and φ 2 which are computed with relation (66) using coefficients CL 1 = 1 and CL 2 = 10. The top diagram presents function φ 1 which is equal to 1 11 between instants 2 and 4 and to 1 between instants 0 and 2 and between instants 4 and 6. The middle diagram presents function φ 2 which is equal to 10 11 between instants 2 and 4 and to 0 elsewhere. The bottom diagram shows that the sum of functions φ 1 and φ 2 is equal to 1 between instants 0 and 6. This example illustrates that formula (66) computes correctly the partition of unity. Now, if the fuzzy index F I is equal to 1 in simulation of Figure 5, functions φ 1 and φ 2 are piecewise rational fractions over time periods [START_REF] Egger | Tikhonov regularization applied to the inverse problem of option pricing: convergence analysis and rates[END_REF][START_REF] Frénod | A continuous-in-time financial model[END_REF] and [START_REF] Frénod | A continuous-in-time financial model[END_REF][START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF] that are in form of fractions of affines functions.

The left picture of Figure 5 shows the partition of unity formed by functions φ 1 and φ 2 which are computed with relation (65). The top diagram presents function φ 1 which is a rational function over time periods [START_REF] Egger | Tikhonov regularization applied to the inverse problem of option pricing: convergence analysis and rates[END_REF][START_REF] Frénod | A continuous-in-time financial model[END_REF] and [START_REF] Frénod | A continuous-in-time financial model[END_REF][START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF] and is equal to 1 2 at instant 3. Function φ 1 is equal to 1 over time periods [0, 2] and [START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF][START_REF] Frńod | Continuous-in-time financial model for public communities[END_REF]. The middle diagram presents function φ 2 which is a rational function over time periods [START_REF] Egger | Tikhonov regularization applied to the inverse problem of option pricing: convergence analysis and rates[END_REF][START_REF] Frénod | A continuous-in-time financial model[END_REF] and [START_REF] Frénod | A continuous-in-time financial model[END_REF][START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF] and is equal to 1 2 at instant 3. Function φ 2 is null over time periods [0, 2] and [START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF][START_REF] Frńod | Continuous-in-time financial model for public communities[END_REF]. The bottom diagram shows that function φ 1 + φ 2 is equal to 1 between instants 0 and 6. Consequently, relation (65) computes correctly the partition of unity.

The right picture of Figure 5 shows the partition of unity formed by functions φ 1 and φ 2 which are computed with relation (66) using coefficients CL 1 = 1 and CL 2 = 10. The top diagram presents

φ1φ2φ1φ2

  The first function φ1 of the partition of unityThe second function φ2 of the partition of unityThe first function φ1 of the partition of unityThe second function φ2 of the partition of unity

Figure 4 :

 4 Figure 4: In each of the two pictures, the fuzzy index F I is equal to 0. In the left picture, function φ 1 is represented at the top, function φ 2 in the middle which are computed with relation (65) and their sum φ 1 + φ 2 at the bottom. In the right picture, function φ 1 is represented at the top, function φ 2 in the middle which are computed with relation (66) with coefficients CL 1 = 1 and CL 2 = 10 and the sum φ 1 + φ 2 at the bottom.

φ1φ2φ1φ2

  The first function φ1 of the partition of unityThe second function φ2 of the partition of unityThe first function φ1 of the partition of unityThe second function φ2 of the partition of unity

Figure 5 :

 5 Figure 5: In each of the two pictures, the fuzzy index F I is equal to 1. In the left picture, function φ 1 is represented at the top, function φ 2 in the middle which are computed with relation (65) and their sum φ 1 + φ 2 at the bottom. In the right picture, function φ 1 is represented at the top, function φ 2 in the middle which are computed with relation (66) with coefficients CL 1 = 1 and CL 2 = 10 and the sum φ 1 + φ 2 at the bottom.

  with FI = 1, CL1 = 1, CL2 = 10

Figure 6 :

 6 Figure 6: In each of the two pictures. In the left picture, fuzzy index F I equals to 0, Loan Measure κE is represented at the top, in the middle the concentrated measure of κE is computed with the partition of unity presented at the left picture of Figure4, in the bottom the concentrated measure of κE is computed with the partition of unity presented at the right picture of Figure4. In the right picture, fuzzy index F I equals to 1, κE is represented at the top, in the middle the concentrated measure of κE is computed with the partition of unity presented at the left picture of Figure5, in the bottom the concentrated measure of κE is computed with the partition of unity presented at the right picture of Figure5.

  If Repayment Pattern Measure γ satisfies relation (11) and if Loan Measure κE is in M([t I , Θ max -Θ γ ]) and Initial Debt Repayment Measure ρI K is in M([t I , Θ max ]), then Algebraic Spending Measure σ is also in M([t I , Θ max ]) and has the next expression in terms of Loan Measure κE :

		has its support in
	[t I , Θ max ], we need to choose Loan Measure κE with support being included in [t I , Θ max -Θ γ ]. Indeed,
	the support of convolution of two compactly supported measures κE and γ is included in the sum of
	their supports which is translated by following relation
	Supp(κ E γ) ⊂ Supp(κ E ) + Supp(γ).	(12)
	Theorem 3.1.

  1 |[tI,t[ dt is density measure which density stands for the function

	t
	t →	κE -κE γ.	(15)
	tI
	Let D be an operator defined in M([t I , Θ max ]) acting on Initial Debt Repayment Measure ρI K by
	D[ρ I K ] = -α ρI K , 1 |[t,Θmax[ dt -ρI K ,	(16)
	where ρI K , 1 |[t,Θmax[ dt is density measure which density stands for the function
		Θmax
	t →	ρI K .	(17)
		t
	Expression (13) of Measure σ can be written as follows:

F I , where point d j is the left neighborhood of the point d i such that a distance d i -d j is strictly positive and is minimized. If there is any point strictly inferior than d i , we take d j = d i . Formally, points (d i ) i∈J are defined as follows:

This work is jointly funded by MGDIS company (http://www.mgdis.fr/) and the LMBA (http://www.lmba-math.fr/).

function φ 1 and the middle diagram presents function φ 2 . Functions φ 1 and φ 2 are rational functions over time periods [START_REF] Egger | Tikhonov regularization applied to the inverse problem of option pricing: convergence analysis and rates[END_REF][START_REF] Frénod | A continuous-in-time financial model[END_REF] and [START_REF] Frénod | A continuous-in-time financial model[END_REF][START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF]. Function φ 1 is equal to 1 over time periods [0, 2] and [START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF][START_REF] Frńod | Continuous-in-time financial model for public communities[END_REF], whereas function φ 2 is null over time periods [0, 2] and [START_REF] Frénod | Optimal control of a continuous-in-time financial model[END_REF][START_REF] Frńod | Continuous-in-time financial model for public communities[END_REF]. The bottom diagram shows that function φ 1 + φ 2 is equal to 1 between instants 0 and 6.

Financial simulation

We will justify, in this section, relation (56), that computes Dirac masses of the concentrated measure. We will adopt both mathematical and financial point of view.

The left picture of Figure 6 represents two concentrated measures of Loan Measure κE . This Loan Measure κE is shown in the top diagram which is a piecewise function that is equal to 1 between instants 0 and 0.5, to 4 between instants 1 and 2, to 6 between instants 3 and 4, to 8 between instants 4 and 5, to 3 between instants 5 and 5.5 and to 0 elsewhere. The first concentrated measure (see middle diagram) is computed with fuzzy index F I equals to 0 . The second concentrated measure (see bottom diagram) is computed with fuzzy index F I equals to 0 and with coefficients CL 1 = 1 and CL 2 = 10. We precise that we use the partition of unity presented at the left picture and at the right picture of Figure 4 to compute respectively the first concentrated measure and the second concentrated measure.

The measure drawn in the middle diagram is computed with formula (56) and is a combination of two Dirac masses (the first Dirac mass has the mass 17 at time 2.5 and the second Dirac mass has the mass 3 at time 3.5). This simulation can be interpreted from the financial slant as follows. It means that the loan is shared into two pieces, the first one consists in borrowing 17 at time 2.5 and the second one consists in borrowing 3 at time 3.5. Thus, the total borrowed amount is equal to 17 + 3 = 20. This makes that formula (58) is homogeneous with regard to this simulation.

The bottom diagram shows the concentrated measure which is computed with relation (56). It is a combination of two Dirac masses (the first Dirac mass has the mass 160 11 at time 2.5 and the second Dirac mass has the mass 60 11 at time 3.5). The meaning of this simulation can be given as follows. It means that the loan is shared into two pieces, the first one consists in borrowing 160 11 at time 2.5 and the second one consists in borrowing 60 11 at time 3.5. This loan is the localized loan of κE at instants 2.5 and 3.5. It means that the total borrowed amount is made of amount 160 11 at times 2.5 and amount 60 11 at time 3.5, which is equal to 160 11 + 60 11 = 20. Consequently, since the amount borrowed associated to Loan Measure κE between times 0 and 6 is 20, relation (58) is homogeneous.

In order to illustrate formulas (56) and (58), in the right picture of Figure 6, we draw examples of concentrated measures of κE with fuzzy index F I equals to 1. This picture shows three measures, the first measure is Loan Measure κE shown in the top diagram which is the same as previously (see top diagram of the left picture). The other measures are the concentrated ones of κE . We use the partition of unity presented at the left picture and at the right picture of Figure 5 to compute respectively the first concentrated measure (see middle diagram) and the second concentrated measure (see bottom diagram).

These two concentrated measures are computed with formula (56). The first concentrated measure in the middle diagram is 17.5δ t=2.5 + 2.5δ t=3. [START_REF] Frénod | Two optimization problems using a continuous-in-time financial model[END_REF] . This illustrates that the loan gives the way that an amount 17.5 is borrowed at time 2.5 and an amount 2.5 is borrowed at time 3.5. It means that the total borrowed amount is equal to 20. The second concentrated measure in the bottom diagram is 15.2δ t=2.5 + 4.8δ t=3. [START_REF] Frénod | Two optimization problems using a continuous-in-time financial model[END_REF] . This concentrated measure means that an amount 15.2 is borrowed at time 2.5 and an amount 4.8 is borrowed at time 3.5. It means that the total borrowed amount is equal to 20. As the amount borrowed associated to Loan Measure κE between times 0 and 6 is equal to 20, formula (58) is correct. It is concluded that this financial simulation using concentrated measures gives realistic results.


Guillaume Favelier
email: guillaume.favelier@lip6.fr

Charles Gueunet
email: charles.gueunet@kitware.com

Julien Tierny
email: julien.tierny@lip6.fr

Visualizing Ensembles of Viscous Fingers

Keywords: Topological data analysis, data segmentation, feature tracking, ensemble data, viscous fingers

This paper presents a topological data analysis framework based on persistent homology and Morse complexes for the visualization and analysis of ensemble data-sets representing viscous fingers. Our approach quantitatively corroborates the classical description of viscous fingers. A systematic analysis across several data resolutions indicates converging statistics as the resolution increases.

INTRODUCTION

Viscous fingering is an instability phenomenon which occurs in porous media at the interface between two fluids of distinct viscosity. In particular, it appears when a less viscous fluid is injected within a more viscous one. It intervenes in many fields of science and engineering [START_REF] De Wit | Viscous fingering of miscible slices[END_REF][START_REF] Feder | Viscous Fingering in Porous Media[END_REF], including geology, hydrology as well as in oil industry where it plays a key role in the extraction process.

This phenomenon is characterized by the formation of a characteristic pattern, called viscous fingers (Fig. 1). In particular, the geometrical evolution of these patterns provides good indications about the evolution of the penetration of the less viscous fluid. Thus, capturing, tracking and analyzing the geometry of viscous fingers is of first importance for the understanding of the penetration process.

Viscous fingering can be decomposed into three major regimes: 1. Launch: Initially, the interface between the two fluids is approximately planar. Then, the less viscous fluid starts to penetrate the more viscous one when sufficient injection force is applied to it. Viscous fingers of low and uniform amplitude start to appear in this phase. 2. Expansion: Once the reaction is launched, the difference of pressure between the two fluids tends to favor an acceleration of the penetration speed for the areas where the less viscous fluid penetrates the most the more viscous fluid. In other words, larger fingers will tend accelerate and grow faster than smaller ones, resulting eventually in a configuration where only large and elongated fingers prevail (Fig. 1). 3. Termination: Optionally, depending on the characteristics of the media and the fluids, the two fluids can eventually mix together in a termination state, where the finger pattern has completely disappeared, after the merge of the large fingers. Depending on the type of data (acquired or simulated), the robust identification and tracking through time of the fingers for each of these steps is subject to several challenges:

1. The represented fluids or the input data representation may not be conducive to an easy discrimination of the less viscous fluid from the other one.

2. Once the fluids have been separated, the shape of their interface can admit a complex geometry, with several protrusions of a variety of scales, which challenges the clear identification of the individual fingers. 3. Once the geometry of the interface has been extracted, the rapid and non-uniform growth of the fingers challenges their robust tracking through time. [START_REF] Floriani | Morse complexes for shape segmentation and homological analysis: discrete models and algorithms[END_REF]. Finally, once a reliable tracking of the fingers has been achieved, relevant quantitative statistics still need to be derived and evaluated for the interpretation of the phenomenon. In this paper, we focus on the viscous fingering simulation data provided in the context of the scientific visualization contest 2016 [START_REF] Ieeevis | Scientific visualization contest[END_REF], which on top of the aforementioned challenges, brings an additional level of complexity to the problem by studying ensemble data. This simulation studies viscous fingering in the context of the mix of continuously dissolving salt and water. As further described in [START_REF] Ieeevis | Scientific visualization contest[END_REF], the simulation runs are given as time-varying particle data-sets representing salt concentration. Due to the stochastic nature of the simulation algorithm, several runs are considered, at distinct resolutions. This yields the following additional challenges:

1. The overall volume of generated data prevents a fast and easy analysis, visualization and interpretation of the phenomenon. 2. The stochastic nature of the simulation code raises the question of the stability of the fingering process across several runs, which needs to be analyzed. 3. The multi-resolution nature of the data raises the question of the convergence of the simulation code, which also needs to be analyzed. In this paper, we address the two sets of challenges described above and we present an interactive framework for the analysis and visualization of ensembles of viscous fingers. In particular, we show how to extend and adapt to viscous fingering the data analysis pipeline proposed by Laney et al. [START_REF] Laney | Understanding the structure of the turbulent mixing layer in hydrodynamic instabilities[END_REF] in the context of the study of the Rayleigh Taylor instability. Finally we report the findings we made by using our data analysis framework.

DATA PRE-PROCESSING

The input data is given by three sets of simulation runs (one set per resolution). Each run is represented by time-varying particles carrying salt concentration. We pre-process each time-step by computing a volumetric interpolation of the particle data onto a 128 3 regular grid. In particular, we use an adaptation of Shepard's kernel method [START_REF] Shepard | A two-dimensional interpolation function for irregularly-spaced data[END_REF] to three-dimensional domains. Next, this regular grid is clipped to the original shape of the domain (a cylinder) to avoid the effects of the boundary on the interpolation (Fig. 2). This results in a tetrahedral mesh that will be the input to our data analysis pipeline (Sec. 3). The advantage of this preprocessing is twofold:

• First, it brings continuity to the data, which enables the usage of a larger class of analysis algorithms. • Second, it allows for truly interactive renderings of each time step based on volume rendering, isosurfaces or color mapped clipped views for instance. This pre-processing addresses task 1 of the contest instructions [START_REF] Ieeevis | Scientific visualization contest[END_REF].

Our data analysis pipeline adapts the approach by Laney et al. [START_REF] Laney | Understanding the structure of the turbulent mixing layer in hydrodynamic instabilities[END_REF] to viscous fingering. It is composed of 5 steps, illustrated in Fig. 3, which are described in the following subsections.

Fluid discrimination (per time-step)

First, for each time step, we separate the dissolving salt from the ambient water. Given the piecewise linear (PL) scalar field f c representing salt concentration and defined on the tetrahedral mesh M representing our input domain, we identify as dissolving salt the sur-level set L + of f c at the isovalue i salt :

L + (i salt ) = {p ∈ M | f c (p) ≥ i salt } (1)
From our experience, we found that an isovalue i salt of 10 gave consistent results along time steps and across runs. Finally, to ignore spurious bubbles, we only consider in the remainder the largest connected component of L + (i salt ), that we note S (Fig. 3(a)).

Finger identification (per time-step)

Next, given the geometry of the dissolving salt S , we aim at extracting in a robust manner the tips of the viscous fingers. To achieve this, we consider as finger tips the points which are locally the furthest away from the top of the domain, where salt is continuously added. Intuitively, this corresponds to salt particles which traveled the furthest from their origin. We model such points as local maxima [START_REF] Banchoff | Critical points and curvature for embedded polyhedral surfaces[END_REF] of the geodesic distance [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF] f d : S → R + from the top of the domain. As shown in Fig. 3(b) (small light green spheres), this strategy identifies as finger tip even slight bumps in the geometry. Thus, we employ persistent homology [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] to filter the maxima. In particular, our experience shows that preserving maxima whose persistence is higher than 10% of the function span provides consistent results along time-steps and across runs (Fig. 3(b), large dark green spheres). In practice, the persistence of the maxima is evaluated by applying the Elder's rule [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] during the expansion (in decreasing order of f d ) of a union-find data structure [START_REF] Cormen | Introduction to Algorithms[END_REF].

Next, given the list of finger tips T previously identified, we compute for each vertex of S , the geodesic distance to its closest finger tip, noted f t : S → R + (Fig. 3(c)). We finally identify as viscous fingers each cell of the Morse complex [START_REF] Floriani | Morse complexes for shape segmentation and homological analysis: discrete models and algorithms[END_REF] of f t (Fig. 3(d)). In practice, such cells are evaluated with a watershed traversal [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF] initiated at each finger tip. Here, like Laney et al. [START_REF] Laney | Understanding the structure of the turbulent mixing layer in hydrodynamic instabilities[END_REF], the usage of the Morse complex enables a complete partitioning of S , with a bijection from the set of partitions to the set of finger tips (in contrast to a Reeb graph [START_REF] Reeb | Sur les points singuliers d'une forme de Pfaff complètement intégrable ou d'une fonction numérique[END_REF] based segmentation). However, in contrast to Laney et al. [START_REF] Laney | Understanding the structure of the turbulent mixing layer in hydrodynamic instabilities[END_REF], the Morse complex is not computed on the interface between the fluids but on the tetrahedral mesh S , which facilitates the evaluation of volumetric indicators. This step partially addresses the task 2 of the contest instructions [START_REF] Ieeevis | Scientific visualization contest[END_REF].

Finger tracking (per run)

Once fingers have been extracted on a per time-step basis, we proceed to their tracking through time with an approach similar to topology based techniques [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF][START_REF] Sohn | Time varying contour topology[END_REF]. In particular, for each finger at a time-step t, we connect it to the finger at time-step t + 1 which maximizes the volume of their intersection (Fig. 3(e)). Once all the time-steps have been processed, we re-iterate through the list of fingers to give a unique identifier to fingers that are connected through time. If multiple fingers of a time step t connect to a single finger at time step t + 1, we give it the identifier of the oldest of the fingers of time step t. This heuristic is motivated by topological persistence [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] (Elder's rule). Finally, we stop the tracking procedure when fingers hit the bottom of the domain. An example of time-tracking is shown in Fig. 4. The output of our time tracking algorithm is a series of tetrahedral meshes S i , with a finger identifier defined per vertex. Each of these can be visualized interactively with our system and the users can interactively navigate through the time steps. This finishes to address task 2 of the contest instructions [START_REF] Ieeevis | Scientific visualization contest[END_REF].

Quantitative analysis (per run)

Once the fingers have been tracked through time, various timevarying statistics are computed on a per finger basis, such as the evolution of its volume for instance (see Sec. 5 for a comprehensive list of measures). Each of these statistics is shown to the users as a 1D plot over time with a color code matching that of the segmentation in the 3D view. Then, if the users identify an irregular pattern in the curve of the time-varying statistics, they can directly inspect the corresponding finger in 3D (Fig. 5) and navigate through the corresponding time steps. We also provide global time-varying statistics computed over the set of fingers, or over the dissolving salt S . As described in Sec. 5, for a given run, users typically first inspect the time-varying 1D plots and then further focus on the corresponding fingers/time steps with the time-varying 3D linked view. This addresses task 3 of the contest instructions [START_REF] Ieeevis | Scientific visualization contest[END_REF].

Comparative analysis (across runs)

Our user interface also offers comparative analysis capabilities by supporting the side by side display of multiple time-varying global statistics (one per run), as well as 3D views that are linked to each of the time-varying statistics window (Fig. 6). Typically, users first inspect a summary view of the time-varying global statistics for all runs (Fig. 7). This enables the visual detection of repeating patterns across runs as well as outliers. Then, our interface enables the interactive inspection in comparative views of the time-varying fingers for the runs selected by the users. This addresses task 4 of the contest instructions [START_REF] Ieeevis | Scientific visualization contest[END_REF].

IMPLEMENTATION

We implemented our data analysis framework in C++. In particular, processing tasks have been designed as VTK filters [START_REF] Kitware | The Visualization Tool Kit (vtk[END_REF] which were globally handled by Python scripts. These scripts were run offline on commodity desktop computers to compute the finger extraction and tracking (Sec. 3.1 through Sec. 3.3) as well as the time-varying statistics (Sec. 3.4 and Sec. 3.5). We wrote ParaView plugins [START_REF]Paraview[END_REF] to implement the user interface of our system.

The typical preprocessing time (Sec. 2) for a single run is approximately 2 hours, while the complete run time for our data analysis framework is approximately one hour per run. We applied our analysis pipeline to the 93 runs provided in the framework of the contest (22 high-res, 23 medium-res and 48 low-res runs). Regarding the exploration features of our framework, these run at interactive rates, as showcased in the companion video, which presents live captures that were shot in real time. This addresses task 5 of the contest instructions [START_REF] Ieeevis | Scientific visualization contest[END_REF].

RESULTS AND DISCUSSION

This section reports the findings we made by using our data analysis framework. In particular, we first focus our analysis on the 22 runs at high resolution (1.7M particles). The comparative analysis across resolutions (544k and 194k particles) is discussed in Sec. 5.3.

Regime identification

In this section, we first try to corroborate the decomposition of the fingering process in three regimes (Sec. 1). To do so, we will first inspect summary views of global statistics (Sec. 3.5) for all runs.

We first analyze the evolution of the descent of the dissolving salt, by looking at the minimum Z-coordinate of S through time (Fig. 7). Three major observations come out of this plot. First, the descent seems linear with time after time-step 10. Second, all runs behave identically in the early time-steps and only little differences can be observed towards the end of the process. Moreover, such a difference seems to increase linearly with time. Third, the 3 regimes identified in Sec. 1 (Launch, Expansion, Termination) are clearly visible: Termination: All runs hit the bottom of the domain in the same time-frame (time step 50 to 60), after which a mixing of the two fluids can be observed (see the companion video). Launch:

From time-step 0 to 10, the dissolving salt remains at the top of the domain; Expansion: From time-step 10 to 50-60, the dissolving salt travels down the domain.

To confirm the bounds of the launch regime, we inspect the evolution of the average salt concentration of the dissolving salt S (Fig. 8). These curve exhibit a clear inflection pattern in the early stages of the simulation, with a clear local maximum at time-step 10. We suspect this maximum marks the end of the launch phase. To confirm this hypothesis, we inspect in comparative mode (Sec. 3.5) the runs where the dissolving salt hits first (Run #19) and last (Run #15) the bottom of the domain. Fig. 9 shows in comparative views the extracted finger tips (green spheres) along with the geodesic distance to the top of the domain ( f d , blue to green color gradient). Note that the geometry of S , the set of extracted finger tips and the level lines of f d are highly similar in both runs in the early stage of the launch regime (time-step 5, top). This similarity quickly vanishes towards the end of the launch regime (time-step 10, bottom), where small amplitude fingers clearly start to appear (small bumps in S , bottom). After this time-step, fingers clearly evolve independently in each run. Thus, we conclude that the launch regimes terminates at time-step 10 for all runs.

Characteristics of the expansion regime

In this section, we try to corroborate the description of the fingering process in the expansion regime (Sec. 1), where larger fingers are supposed to grow faster, at the expense of smaller ones, which they eventually absorb. To do so, we first inspect the evolution of the number of fingers through time (Fig. 10). This plot confirms that the number of fingers globally decreases for all runs, with a consistent decrease rate across runs, to eventually tend to a small number (typically five or less) towards the end of the expansion regime. Fig. 11 provides a comparative view for the fastest and slowest runs, which confirms the existence of only few, very large fingers at the end of the expansion regime. To further understand how fingers evolve during the expansion process, we inspect the evolution through time of the ratio between the volume of the largest finger and that of the dissolving salt S (Fig. 12). These curve indicate how the largest finger evolves in size for all runs. In particular, Runs #19 and #4 show a particularly large finger towards the end of the expansion process. We will focus on these to better understand the decrease in the number of fingers. Figs. 13 and14 show comparative views before and after the merge of large fingers in Run #19 and #4 respectively. In particular, these figures show that the fingers that are further down the domain are also the fastest.

By further inspecting other finger merge events, we identified an unexpected behavior, that we call detachment, (Figs. 15 and16 for Runs #19 and #4 respectively). In particular, these figures show that clearly prominent fingers, which are both the fastest and the deepest down the domain, can suddenly detach from the rest of the dissolving salt S . This phenomenon is preceded by a sudden gain of velocity, which is accompanied by a dramatic drop down in salt concentration, eventually fully dissolving the finger in the ambient water. Our inspections reveal that such detachments occur several times in all runs, either early or late in the expansion. We suspect these may have an impact on the efficiency of the penetration of the less viscous fluid since the curves of average salt concentration indicate that detaching fingers instantly vanish and dissolve in the ambient water, suddenly stopping the progression of the dissolving salt in these regions (see Sec. 6 for further comments).

Inter-resolution analysis

In this section, we study the impact of the input data particle resolution on the characterization of the viscous fingering process. To do so, we visualize global summarizes provided by global time-varying statistics for all runs in Fig. 17. This figure plots, for the dissolving salt S , the evolution over time of the minimimum Z-coordinate, the average salt concentration, the volume, the average velocity and the number of fingers. In particular, each of the three resolutions is represented with a distinct color. These three sets of curves exhibit similar global behaviors: a salt descent according to a linear slope (Fig. 17(a)), a linear decrease of salt concentration in the expansion regime (Fig. 17(b)) and a linear increase in both volume (Fig. 17(c)) and velocity (Fig. 17(d)). However, each resolution can be easily distinguished from the others as curves of the same color tend to cluster, with only few overlap with the other colors. This indicates clear distinctions between resolutions (see Sec. 6 for further comments).

CONCLUSION

In this paper, we have presented a topological data analysis framework, inspired by Laney et al. [START_REF] Laney | Understanding the structure of the turbulent mixing layer in hydrodynamic instabilities[END_REF], that we adapted to the special case of the analysis of ensembles of particle data-sets representing viscous fingers. Our analysis pipeline is based on a volumetric interpolation of the particle data (Sec. 2), followed by a finger tip extraction based on persistent homology and a finger identification with the Morse complex (Sec. 3.2). The tracking through time of the fingers by overlap estimation (Sec. 3.3) enables the evaluation of relevant time-varying statistics (Secs. 3.4 and 3.5).

By using our data analysis framework on the data-sets provided in the context of the scientific visualization contest 2016 [START_REF] Ieeevis | Scientific visualization contest[END_REF], we were able to corroborate the characterization of the viscous fingering process into three distinct regimes (launch, expansion, termination) for which we were able to specify the exact intervals of appearance. The interactive exploration features of our framework also enabled the quantitative validation of the characterization of the expansion regime, where the deepest fingers indeed grow the fastest, at the expense of the smaller ones that they eventually absorb. This explains the consistent decrease of the number of fingers through time. Additionally, we observed an unexpected behavior, that we called detachment, where fast, deep and with low salt concentration fingers suddenly detach from the rest of the dissolving salt. Our inter-resolution analysis revealed that simulation runs of distinct resolutions had globally similar behavior. However, later expansion starts, faster penetration rates and fewer viscous fingers were observed as the resolution decreases. This comparison shows that the similarity to the highest resolution runs (according to the above criteria) decreases with the resolution, suggesting a convergence of the simulation code for increasing resolutions.

In the future, we would like to validate our analysis framework with domain experts. Moreover, we would like to investigate the design of quantitative measures for a deeper evaluation of detachments, as we suspect they may have an impact on the efficiency of the penetration. The 3D linked views show a clear prominence of the dark blue finger, despite a the very thin portion that connects it to the rest of the dissolving salt S . This portion becomes so thin at time-step 32 that this finger gets disconnected from S and is therefore no longer considered as a finger by our analysis pipeline. We call such a phenomenon detachment, since a rapid finger grows so fast that it detaches from the rest of dissolving salt S and rapidly vanishes in the ambient water. Time varying statistics are provided on the right side: minimum of Z-coordinate for each finger, average velocity for each finger and average salt concentration for each finger. These curves are shown with dash lines for fingers which are no longer active at the current time step. The bold black dash lines indicate global statistics for the region of dissolving salt S . The bold, dark blue curves indicate time-varying statistics for the dark blue finger. In particular, these curve show that, before detachment, this finger was both the deepest down the domain (top curves) and the fastest (center curves). Moreover, the bottom curves show that this finger had a dramatic drop down of salt concentration (bottom curves) exactly when it started to accumulate speed (time-step 27), to eventually fully dissolve in the ambient water.

Figure 16: Finger detachment (dark red finger) at time-step 43 of Run #4 (top: before, bottom: after). The 3D linked views show a clear prominence of the dark red finger, despite a the very thin portion that connects it to the rest of the dissolving salt S . This portion becomes so thin at time-step 43 that this finger gets disconnected from S and is therefore no longer considered as a finger by our analysis pipeline. We call such a phenomenon detachment, since a rapid finger grows so fast that it detaches from the rest of dissolving salt S and rapidly vanishes in the ambient water. Time varying statistics are provided on the right side: minimum of Z-coordinate for each finger, average velocity for each finger and average salt concentration for each finger. These curves are shown with dash lines for fingers which are no longer active at the current time step. The bold black dash lines indicate global statistics for the region of dissolving salt S . The bold, dark red curves indicate time-varying statistics for the dark red finger. In particular, these curve show that, before detachment, this finger was both the deepest down the domain (top curves) and the fastest (center curves). Moreover, the bottom curves show that this finger had a dramatic drop down of salt concentration (bottom curves) exactly when it started to accumulate speed (time-step 38), to eventually fully dissolve in the ambient water. (e) The evolution of the number of fingers indicates a clear distinction for the lowest resolution in the early time-steps. However, this increase is not confirmed in the expansion regime (typically around time-step 25), where fewer fingers are extracted as the resolution decreases. In conclusion, runs of distinct resolutions exhibit similar global behavior, however with later expansion starts, faster penetration rates and fewer viscous fingers as the resolution decreases. This comparative study shows that the similarity to the highest resolution runs (according to the above criteria) decreases with the resolution, suggesting a convergence of the simulation code for increasing resolutions.

Figure 1 :

 1 Figure 1: Photograph of viscous fingers in a glycerol-water mix (creative commons license).

Figure 2 :Figure 3 :

 23 Figure 2: Preprocessing for each time-step. The input particle data is shown with small spheres (screen space rendering) and the salt concentration is shown with the color gradient (a). Particles are interpolated onto a regular grid with Shepard's method [14] (b). Next, the grid is clipped (c) to the original domain's shape (a cylinder) to avoid boundary effects. (d) and (e) show zoom-ins for the white square ((a and (b)). Data-sets are shown upside down to reduce occlusion.

Figure 4 :

 4 Figure 4: Example of finger tracking through time on the run #3 of the lowest particle resolution. Fingers that are connected through time are shown with the same color.

Figure 5 :

 5 Figure5: Screen-shot of the user interface to our data-analysis framework (intra-run mode). Users can select from the bottom-left list the per-finger, time-varying statistics to visualize on the right side of the screen (white background). There, the vertical black line indicates the time-step being currently visualized in the linked 3D view (center). Users can navigate through time steps with the time navigation buttons (top, center). Various per-finger statistics can be displayed in the 3D view by pointing on a finger with the cursor (dark rectangle).

Figure 6 :

 6 Figure 6: Screen-shot of the user interface to our data-analysis framework (comparative inter-run mode). Users can select from the bottom-left list the global time-varying statistics to visualize on the right side of the screen. In this case, the evolution of the number of fingers, of the fingers' volume, of average salt concentration and average velocity are shown for the runs #1 and #3 of the lowest particle resolution. There, the vertical black line indicates the time-step being currently visualized in the linked 3D views (center top: run #3, center bottom: run #1). Users can navigate through time steps with the time navigation buttons (top, center).

Figure 7 :

 7 Figure 7: Evolution of the minimum Z-coordinate of the region of dissolving salt S for the 22 runs at high resolution. The descent evolution is linear through time (from time-step 10 to time-steps 50-60) and it is consistent across runs. For all runs, the dissolving salt is hitting the bottom of the domain in between time-steps 50 and 60.

Figure 8 :

 8 Figure 8: Evolution of the average salt concentration through time for the region of dissolving salt S , for the 22 runs at high resolution. These curves exhibit a clear inflection pattern in the early stages of the simulation, with a clear local maximum of average concentration at time-step 10. This maximum seems to mark the end of the launch phase.

Figure 9 :

 9 Figure 9: Interactive inspection of the launch regime in early and late time-steps (top: time-step 5, bottom: time-step 10) for the fastest and slowest runs (Run #19 and #15 respectively). The time-varying statistics window (right) plots the evolution of salt concentration in the region of dissolving salt S . Despite the important difference in termination times (time-step 51 vs 60), both runs exhibit very similar patterns in the launch regime. The left 3D linked views show the corresponding regions S , with the extracted finger tips (green spheres) and the geodesic distance to the top of the domain ( f d , blue to green color gradient). The geometry of S , the set of extracted finger tips and the level lines of f d look alike in both runs in the early stages of the launch regime (time-step 5, top). However, this similarity quickly vanishes towards the end of the launch regime (time-step 10, bottom), where small-amplitude fingers clearly start to appear (see the small bumps in S in the linked 3D views).

Figure 10 :

 10 Figure10: Evolution of the number of fingers through time for the 22 runs of the highest resolution. For all runs, this number globally decreases and converges towards a low number (typically five or less) towards the end of the expansion regime. Also note that the decrease rate is consistent across runs (similar slopes).

Figure 11 :

 11 Figure 11: Extracted fingers at the end of the expansion regime for the fastest and slowest runs (Run #19 and #15 respectively) of the highest resolution. The 3D linked views indicate large prominent fingers (one color per finger) in both runs.

Figure 12 :

 12 Figure 12: Evolution through time of the ratio between the volume of the largest finger and that of the dissolving salt S , for the 22 runs of the highest resolution. These curve indicate how the largest finger evolves in size for all runs. In particular, Runs #19 and #4 show a particularly large finger towards the end of the expansion (time-steps 50, purple curve, and 54, dark green curve, respectively).

Figure 13 :

 13 Figure13: Merging of two large fingers in Run #19, before (top) and after (bottom). The 3D linked views show that the prominence of the red finger is not sufficient anymore for it to be identified as an independent finger after time-step 48. Technically, this is due to the fact that the persistence of the corresponding finger tip went below the admissible threshold. At this point, the cyan and red fingers are considered as merging. Time varying statistics are provided on the right side: minimum of Z-coordinate for each finger, average velocity for each finger and average salt concentration for each finger. These curves are shown with dash lines for fingers which are no longer active at the current time step. The bold black dash lines indicate global statistics for the region of dissolving salt S . The bold, red curves indicate time-varying statistics for the red finger. In particular, these curve show that, before the merging, this finger was both the deepest down the domain (top curves) and the fastest (center curves).

Figure 14 :

 14 Figure14: Merging of two large fingers in Run #4, before (top) and after (bottom). The 3D linked views show that the prominence of the red finger is not sufficient anymore for it to be identified as an independent finger after time-step 48. Technically, this is due to the fact that the persistence of the corresponding finger tip went below the admissible threshold. At this point, the yellow and red fingers are considered as merging. Time varying statistics are provided on the right side: minimum of Z-coordinate for each finger, average velocity for each finger and average salt concentration for each finger. These curves are shown with dash lines for fingers which are no longer active at the current time step. The bold black dash lines indicate global statistics for the region of dissolving salt S . The bold, yellow curves indicate time-varying statistics for the yellow finger. In particular, these curve show that, before the merging, this finger was both the deepest down the domain (top curves) and the fastest (center curves).

Figure 15 :

 15 Figure15: Finger detachment (dark blue finger) at time-step 32 of Run #19 (top: before, bottom: after). The 3D linked views show a clear prominence of the dark blue finger, despite a the very thin portion that connects it to the rest of the dissolving salt S . This portion becomes so thin at time-step 32 that this finger gets disconnected from S and is therefore no longer considered as a finger by our analysis pipeline. We call such a phenomenon detachment, since a rapid finger grows so fast that it detaches from the rest of dissolving salt S and rapidly vanishes in the ambient water. Time varying statistics are provided on the right side: minimum of Z-coordinate for each finger, average velocity for each finger and average salt concentration for each finger. These curves are shown with dash lines for fingers which are no longer active at the current time step. The bold black dash lines indicate global statistics for the region of dissolving salt S . The bold, dark blue curves indicate time-varying statistics for the dark blue finger. In particular, these curve show that, before detachment, this finger was both the deepest down the domain (top curves) and the fastest (center curves). Moreover, the bottom curves show that this finger had a dramatic drop down of salt concentration (bottom curves) exactly when it started to accumulate speed (time-step 27), to eventually fully dissolve in the ambient water.

Figure 17 :

 17 Figure 17: Global statistics as a function of time for the 22 high-res (green curves), 23 medium-res (red curves) and 48 low-res (blue curves) runs. These three sets of curves exhibit similar global behaviors: a salt descent according to a linear slope (Fig. 17(a)), a linear decrease of salt concentration in the expansion regime (Fig. 17(b)) and a linear increase in both volume (Fig.17(c)) and velocity (Fig. 17(d)). However, each resolution can be easily distinguished from the others as curves of the same color tend to cluster, with only few overlap with the other colors. This indicates clear distinctions between resolutions. (a) The minimum Z-coordinate of the dissolving salt S indicates clearly distinct descent rate for the deepest fingers (different slopes). Here lower resolutions hit the bottom of the domain faster. (b) The average salt concentration of the dissolving salt S indicates clearly distinct initial salt concentration levels, with a slight delay for the start of the expansion regime (local maximum of concentration) as the resolution decreases. Note that from the middle of the expansion regime onward, the decrease rate is similar across resolutions. (c) The volume of the dissolving salt S also indicates different growth rate (lower resolutions grow faster, as suggested by (a)). (d) The average velocity within the dissolving salt S also indicates different average speeds, confirming the increase in speed for lower resolutions.(e) The evolution of the number of fingers indicates a clear distinction for the lowest resolution in the early time-steps. However, this increase is not confirmed in the expansion regime (typically around time-step 25), where fewer fingers are extracted as the resolution decreases. In conclusion, runs of distinct resolutions exhibit similar global behavior, however with later expansion starts, faster penetration rates and fewer viscous fingers as the resolution decreases. This comparative study shows that the similarity to the highest resolution runs (according to the above criteria) decreases with the resolution, suggesting a convergence of the simulation code for increasing resolutions.

  Figure 17: Global statistics as a function of time for the 22 high-res (green curves), 23 medium-res (red curves) and 48 low-res (blue curves) runs. These three sets of curves exhibit similar global behaviors: a salt descent according to a linear slope (Fig. 17(a)), a linear decrease of salt concentration in the expansion regime (Fig. 17(b)) and a linear increase in both volume (Fig.17(c)) and velocity (Fig. 17(d)). However, each resolution can be easily distinguished from the others as curves of the same color tend to cluster, with only few overlap with the other colors. This indicates clear distinctions between resolutions. (a) The minimum Z-coordinate of the dissolving salt S indicates clearly distinct descent rate for the deepest fingers (different slopes). Here lower resolutions hit the bottom of the domain faster. (b) The average salt concentration of the dissolving salt S indicates clearly distinct initial salt concentration levels, with a slight delay for the start of the expansion regime (local maximum of concentration) as the resolution decreases. Note that from the middle of the expansion regime onward, the decrease rate is similar across resolutions. (c) The volume of the dissolving salt S also indicates different growth rate (lower resolutions grow faster, as suggested by (a)). (d) The average velocity within the dissolving salt S also indicates different average speeds, confirming the increase in speed for lower resolutions.(e) The evolution of the number of fingers indicates a clear distinction for the lowest resolution in the early time-steps. However, this increase is not confirmed in the expansion regime (typically around time-step 25), where fewer fingers are extracted as the resolution decreases. In conclusion, runs of distinct resolutions exhibit similar global behavior, however with later expansion starts, faster penetration rates and fewer viscous fingers as the resolution decreases. This comparative study shows that the similarity to the highest resolution runs (according to the above criteria) decreases with the resolution, suggesting a convergence of the simulation code for increasing resolutions.

  Figure 17: Global statistics as a function of time for the 22 high-res (green curves), 23 medium-res (red curves) and 48 low-res (blue curves) runs. These three sets of curves exhibit similar global behaviors: a salt descent according to a linear slope (Fig. 17(a)), a linear decrease of salt concentration in the expansion regime (Fig. 17(b)) and a linear increase in both volume (Fig.17(c)) and velocity (Fig. 17(d)). However, each resolution can be easily distinguished from the others as curves of the same color tend to cluster, with only few overlap with the other colors. This indicates clear distinctions between resolutions. (a) The minimum Z-coordinate of the dissolving salt S indicates clearly distinct descent rate for the deepest fingers (different slopes). Here lower resolutions hit the bottom of the domain faster. (b) The average salt concentration of the dissolving salt S indicates clearly distinct initial salt concentration levels, with a slight delay for the start of the expansion regime (local maximum of concentration) as the resolution decreases. Note that from the middle of the expansion regime onward, the decrease rate is similar across resolutions. (c) The volume of the dissolving salt S also indicates different growth rate (lower resolutions grow faster, as suggested by (a)). (d) The average velocity within the dissolving salt S also indicates different average speeds, confirming the increase in speed for lower resolutions.(e) The evolution of the number of fingers indicates a clear distinction for the lowest resolution in the early time-steps. However, this increase is not confirmed in the expansion regime (typically around time-step 25), where fewer fingers are extracted as the resolution decreases. In conclusion, runs of distinct resolutions exhibit similar global behavior, however with later expansion starts, faster penetration rates and fewer viscous fingers as the resolution decreases. This comparative study shows that the similarity to the highest resolution runs (according to the above criteria) decreases with the resolution, suggesting a convergence of the simulation code for increasing resolutions.

  Figure 17: Global statistics as a function of time for the 22 high-res (green curves), 23 medium-res (red curves) and 48 low-res (blue curves) runs. These three sets of curves exhibit similar global behaviors: a salt descent according to a linear slope (Fig. 17(a)), a linear decrease of salt concentration in the expansion regime (Fig. 17(b)) and a linear increase in both volume (Fig.17(c)) and velocity (Fig. 17(d)). However, each resolution can be easily distinguished from the others as curves of the same color tend to cluster, with only few overlap with the other colors. This indicates clear distinctions between resolutions. (a) The minimum Z-coordinate of the dissolving salt S indicates clearly distinct descent rate for the deepest fingers (different slopes). Here lower resolutions hit the bottom of the domain faster. (b) The average salt concentration of the dissolving salt S indicates clearly distinct initial salt concentration levels, with a slight delay for the start of the expansion regime (local maximum of concentration) as the resolution decreases. Note that from the middle of the expansion regime onward, the decrease rate is similar across resolutions. (c) The volume of the dissolving salt S also indicates different growth rate (lower resolutions grow faster, as suggested by (a)). (d) The average velocity within the dissolving salt S also indicates different average speeds, confirming the increase in speed for lower resolutions.(e) The evolution of the number of fingers indicates a clear distinction for the lowest resolution in the early time-steps. However, this increase is not confirmed in the expansion regime (typically around time-step 25), where fewer fingers are extracted as the resolution decreases. In conclusion, runs of distinct resolutions exhibit similar global behavior, however with later expansion starts, faster penetration rates and fewer viscous fingers as the resolution decreases. This comparative study shows that the similarity to the highest resolution runs (according to the above criteria) decreases with the resolution, suggesting a convergence of the simulation code for increasing resolutions.

  Figure 17: Global statistics as a function of time for the 22 high-res (green curves), 23 medium-res (red curves) and 48 low-res (blue curves) runs. These three sets of curves exhibit similar global behaviors: a salt descent according to a linear slope (Fig. 17(a)), a linear decrease of salt concentration in the expansion regime (Fig. 17(b)) and a linear increase in both volume (Fig.17(c)) and velocity (Fig. 17(d)). However, each resolution can be easily distinguished from the others as curves of the same color tend to cluster, with only few overlap with the other colors. This indicates clear distinctions between resolutions. (a) The minimum Z-coordinate of the dissolving salt S indicates clearly distinct descent rate for the deepest fingers (different slopes). Here lower resolutions hit the bottom of the domain faster. (b) The average salt concentration of the dissolving salt S indicates clearly distinct initial salt concentration levels, with a slight delay for the start of the expansion regime (local maximum of concentration) as the resolution decreases. Note that from the middle of the expansion regime onward, the decrease rate is similar across resolutions. (c) The volume of the dissolving salt S also indicates different growth rate (lower resolutions grow faster, as suggested by (a)). (d) The average velocity within the dissolving salt S also indicates different average speeds, confirming the increase in speed for lower resolutions.(e) The evolution of the number of fingers indicates a clear distinction for the lowest resolution in the early time-steps. However, this increase is not confirmed in the expansion regime (typically around time-step 25), where fewer fingers are extracted as the resolution decreases. In conclusion, runs of distinct resolutions exhibit similar global behavior, however with later expansion starts, faster penetration rates and fewer viscous fingers as the resolution decreases. This comparative study shows that the similarity to the highest resolution runs (according to the above criteria) decreases with the resolution, suggesting a convergence of the simulation code for increasing resolutions.

ACKNOWLEDGMENTS

This work is partially supported by the Bpifrance grant "AVIDO" (Programme d'Investissements d'Avenir, reference P112017-2661376/DOS0021427) and by the French National Association for Research and Technology (ANRT), in the framework of the LIP6 -Kitware SAS CIFRE partnership reference 2015/1039.