Jean-Noël Jaubert

Laurent Avaullee

Christophe Pierre

Is It Still Necessary to Measure the Minimum Miscibility Pressure?

Gas injection processes are among the most effective methods for enhanced oil recovery. A key parameter in the design of a gas injection project is the minimum miscibility pressure (MMP), the pressure at which the local displacement efficiency approaches 100%. From an experimental point of view, the MMP is routinely determined by slim tube displacements. However, because such experiments are very expensive (timeconsuming), the question we want to answer in this article is as follows: Is this still necessary to measure the MMP? In other words, may other quicker, easier and cheaper gas injection experiments such as swelling test or multicontact test (MCT) substitute for slim tube test? This paper concludes that when the injected gas is not pure CO 2 (and probably not pure N 2 or pure H 2 S), it is enough to fit only two parameters of the equation of state on data including classical PVT data + swelling data + MCT data and then to predict the MMP. The accuracy obtained is similar to the experimental uncertainty. It is thus possible to conclude that the slim tube test may be replaced by swelling tests and MCT, which are much cheaper.

Introduction

Gas injection processes are among the most effective methods for enhanced oil recovery. A key parameter in the design of a gas injection project is the minimum miscibility pressure (MMP), the pressure at which the local displacement efficiency approaches 100%. If the flow is one-dimensional, and if there is no dispersive mixing, then the displacement efficiency is exactly 100% at the MMP. From an experimental point of view, the MMP is routinely determined by slim tube displacements. In those experiments a long (say 10 m or longer), small diameter (say 0.5 cm) tube packed with sand or glass beads is filled with oil that is then displaced by injection gas at a fixed temperature and pressure. The fraction of the initial oil in place recovered after injection of some fixed amount of gas (usually 1.1 or 1.2 pore volume) is measured. The fraction of oil recovered is then plotted as a function of pressure. Typically, recovery increases rapidly with increasing pressure and then levels off. The MMP is usually taken to be the intersection of lines drawn through recovery points in the steeply climbing and level regions as long as the recovery in the level region is above some arbitrary cutoff (often 90%). To determine accurately the MMP, it is in general necessary to perform six displacements at six different pressures. The time required to perform a displacement is about8h(aworking day). This means 1 week of experimental work is necessary to measure one MMP. The cost of such a work depends on the salary of the engineer who makes the experiment but may be evaluated to 10000 US dollars. It is thus clear that for petroleum companies the slim tube test is a very expensive experiment.

The question we want to answer in this article is as follows: Is it still necessary to measure the MMP? To answer yes to this question, two tracks were explored in this paper. The first possibility consists of developing a purely predictive characterization of the heavy cuts allowing an accurate estimation of the MMP. This is why in the first part of this article a purely predictive C 11 + characterization, which we recently developed, is tested. This characterization is perfectly able to predict with a good accuracy all the classical PVT experiments. However, because the MMP is a critical pressure and because the cubic equation of state (EOS) often fails to predict critical points of complex mixtures, it is not always possible to predict accurately the MMP. The deviation between the predicted and the experimental MMP varies between 2% and 20%.

When the first possibility fails, another solution is to fit some parameters of the EOS by using experiments which are much less expensive than the slim tube test (e.g., relative volumes during a constant mass expansion, reservoir fluid bubble pressure, relative volumes during a differential vaporization, stock tank oil density, swelling test data, multicontact test data) and then to predict the MMP with the obtained parameters. The key questions we need to answer are (1) which experiments influence the most the value of the MMP and (2) is it possible to predict the MMP after fitting some parameters on these key experiments. This paper will show that a good fit on the swelling test and above all on the multicontact test (MCT) allows a good prediction of the MMP.

By the end, because the swelling test and MCT were not always conducted experimentally, software able to include the MMP in the fitting procedure was developed. Our PVT package is, to our knowledge, the first one in the world that is able to perform such calculations.

The C 11 + Characterization Procedure

The compositional model used in this paper is a slightly modified version of a previously published work. [START_REF] Avaulle ´e | A New Approach in Correlating the Thermodynamic Oil Properties[END_REF] Indeed, in comparison to our previous work, [START_REF] Avaulle ´e | A New Approach in Correlating the Thermodynamic Oil Properties[END_REF] it was decided to lump together the nine cuts C k (k ) 11-19) in a unique pseudocomponent called "C 11 -C 19 "i n order to reduce the number of components. Indeed, we noticed that such a lumping did not change the predictive power of our compositional model. To do so, the following formulas (well known in petroleum companies), in which x i is the internal mole fraction of cut C i in the pseudocomponent C 11 -C 19 , are used: Moreover, knowing the amount of each cut C k (k ) 11-19) in the reservoir fluid and the groups appearing in each cut C k , it is possible to compute the groups appearing in the pseudocomponent C 11 -C 19 . These groups are necessary to estimate the binary interaction parameters. 1

Predictive Results on 13 Different Crude Oils

TOTALFINAELF performed 15 slim tube tests on 13 reservoir crude oils coming from different countries of the world (we are not allowed to give the name of the countries). For two crude oils, two different gases were injected. The compositional model described in the previous section was used to predict not only these 15 MMPs but also many other measured PVT data. It is here very important to outline that the proposed compositional model is purely predictive in the sense that the only data necessary to perform a flash calculation are, in addition to the reservoir fluid composition, the experimental molar weight and density of each cut from C 11 to C 20 +. No PVT data are used to estimate the parameters of the EOS. Such a model is able to predict:

1. The reservoir fluid bubble pressure with an average deviation of 3.7% (19 bubble point pressures were measured). For most of the crude oils, the bubble point is only measured at the reservoir temperature, but for two of them, such a pressure was measured at three other temperatures.

2. The stock tank oil density with an average deviation of 1.5% (12 values were determined).

3. The evolution of the relative volume during an isothermal constant mass expansion with an average deviation of 0.8% (197 data points were measured). Let us here recall that the relative volume is defined as V rel ) V total /V sat , where V total is the total volume occupied by the fluid at pressure P and V sat is the liquid volume at the bubble pressure.

4. The evolution of the relative volume during an isothermal differential vaporization with an average deviation of 1.5% (71 data points were measured). In this case, the relative volume is defined by V rel ) V liq /V sat , where V liq is the liquid-phase volume at pressure P and V sat is the liquid volume at the bubble pressure.

5. The evolution of the saturation pressure versus quantity of injected gas during a swelling test with an average deviation of 5.4% (42 data points were measured).

6. The evolution of the swollen volume versus quantity of injected gas during a swelling test with an average deviation of 0.7% (42 data points were measured). The swollen volume is defined by V swell ) V sat (k)/V sat0 , where V sat (k) is the volume of the saturated oil at P ) P sat (k), i.e., after the injection of k predetermined amounts of gas. V sat0 is the saturated crude oil volume (before dissolution of any quantity of gas).

7. The evolution of the density of the saturated swollen crude oil (F sat ) versus quantity of injected gas with an average deviation of 1.5%. Forty two data points were measured. More detailed information concerning the swelling test may be found elsewhere. [START_REF] Jaubert | Characterization of Heavy Oils. 3. Prediction of Gas Injection Behavior: Swelling Test, Multicontact Test, Multiple Contact Minimum Miscibility Pressure and Multiple Contact Minimum Miscibility Enrichment[END_REF] 8. The evolution of the partial volume during a MCT (multicontact test) with an average deviation of 24% (13 data points were experimentally determined). In a MCT experiment, a predetermined amount of crude oil and of injection gas are mixed in order to obtain a two-phase system (see ref 2 for more details). The partial volume (PVol) is the ratio of the liquid phase volume by the total cell volume. Although the partial volumes are poorly predicted, our compositional model is able to predict with high accuracy the composition of the phases (liquid and gas) in equilibrium (these results are not shown). 9. The MMP with an average deviation of 7.6% (15 slim tube tests were performed). The algorithm used to calculate the MMP and which takes into account the mixed condensing/vaporizing (C/V) mechanism discovered by Zick [START_REF] Zick | A Combined Condensing/Vaporizing Mechanism in the Displacement of Oil by Enriched Gases[END_REF] is the one developed by Jaubert et al. [START_REF] Jaubert | A Very Simple Multiple Mixing Cell Calculation to Compute the Minimum Miscibility Pressure Whatever the Displacement Mechanism[END_REF] Though our compositional model describes a crude oil with 29 components, we can imagine that the last version of the algorithm developed by Orr and coworkers [START_REF] Wang | Calculation of minimum miscibility pressure[END_REF] could also be used. A discussion on the different algorithms allowing computation the MMP may be found in the articles by Jaubert et al. [START_REF] Jaubert | A Very Simple Multiple Mixing Cell Calculation to Compute the Minimum Miscibility Pressure Whatever the Displacement Mechanism[END_REF][START_REF] Jaubert | Properly Defining the Classical Vaporizing and Condensing Mechanisms When a Gas is Injected Into a Crude Oil[END_REF] Briefly, our algorithm [START_REF] Jaubert | A Very Simple Multiple Mixing Cell Calculation to Compute the Minimum Miscibility Pressure Whatever the Displacement Mechanism[END_REF] is based on a cell to cell simulation with very simple flow dynamics (a P/T flash algorithm is only required to compute the MMP). The recovery factors are calculated at three different pressures correctly selected using 500 cells. To eliminate the numerical dispersion, the results are extrapolated to an infinite number of cells. Typically, the recovery factors increase exponentially with pressure. In our algorithm, the MMP is defined as the pressure at which the recovery factor is equal to 97%.

Table 1 gives a summary of all these results. The numerical values of most of the experimental data used in this work and necessary for a researcher in order to repeat the calculations or compare the results with other methods are published elsewhere. [START_REF] Jaubert | A Crude Oil Data Bank Containing More Than 5000 PVT and Gas Injection Data[END_REF] Figure 1 illustrates graphically the results obtained. The compositions of the 13 reservoir crude oils and of the 17 injected gases are given in Table 2. The gas number, on the first line in Table 2 refers to the crude oil number in which it was injected. As an example, gases G 6a and G 6b were injected in reservoir fluid F 6 . Gases G 1 and G 2 and G 7 and G 8 are identical. It was however decided to give two different numbers to these gases in order to identify simply the crude oil in which they were injected. Gases 5b, 5c, and 5d are quasi-identical. The first one (G 5b ) was used for a swelling test, the second one (G 5c ) for a multicontact test (MCT), and the last one (G 5d ) for a slim tube test. These three gases, which were injected in reservoir fluid F 5 , should be identical. However, because they were prepared at different times, there exists a slight difference between them. This is absolutely not a problem.

The last line in Table 1 shows evidence that the model described in this paper is able to predict the phase behavior of real crude oils with high accuracy. In our previous paper, [START_REF] Avaulle ´e | A New Approach in Correlating the Thermodynamic Oil Properties[END_REF] it is shown that the results obtained with the proposed compositional model compare competitively with those given by other characterizations currently used in petroleum simulators. However, the MMP is predicted with an average deviation of 7.6%, which is twice higher than the experimental uncertainty. To reduce this deviation, a fitting strategy is proposed.

Fitting Strategy To Reduce the Deviations between Experimental and Calculated MMP

Table 1 shows evidence that 6 MMP values are predicted with a deviation lower than 4%, which is the experimental uncertainty (F1 + G1, F2 + G2, F3 + G3, F4 + G4, F5 + G5a, F6 + G6a). For these six cases, it is thus not necessary to proceed to any fitting procedure since the predictive model leads to very nice results. It is here interesting to outline that for the systems F3 + G3 and F5 + G5a, the deviations on the MMP are very small but the deviations on the swelling data are also very small. We may thus suspect a correlation between the MMP and the swelling data (see below).

Concerning the nine other systems, they have to be divided in two families.

The first family concerns the fluids for which the MMP was the unique injection gas experiment performed (no swelling data and no MCT data are available). The five systems concerned are F7 + G7, F8 + G8, F9 + G9, F10 + G10, and F11 + G11. For these five systems, all the available PVT data (relative volumes, reservoir bubble pressure, stock tank oil density) are accurately predicted. It is thus obvious that by fitting the EOS parameters on these experimental data would not lead to a better prediction of the MMP. Indeed the fitted parameters would be very close to the initial ones (given by the predictive procedure). The unique thing to do is thus to fit the parameters of the EOS on all the experimental data, i.e., to include the MMP in the fitting procedure. This is not an easy task because the MMP is calculated using a simplified slim tube in which numerical dispersion is eliminated by extrapolating the results to an infinite number of cells. [START_REF] Jaubert | A Very Simple Multiple Mixing Cell Calculation to Compute the Minimum Miscibility Pressure Whatever the Displacement Mechanism[END_REF] To our knowledge, our software is the first one in the world able to fit the parameters of the EOS (T c , P c , ω, bips) on the MMP and to simultaneously take into account the mixed condensing/vaporizing mechanism. In this study, only two parameters were fitted: the critical temperature and the critical pressure of the pseudocomponent C 11 -C 19 . Such a component was selected because its physical parameters are not wellknown and because its amount in the reservoir fluid is quite high. Table 3 shows evidence that the results obtained are extremely accurate. Indeed this fitting procedure made it possible to estimate the MMP with an average deviation of 2.0%. At the same time, the other PVT properties are also reproduced accurately.

The second family concerns systems for which swelling tests and/or MCT experiments were performed in addition to the slim tube test. The four concerned systems are F5 + G5b, F6 + G6b, F12 + G12, and F13 + G13. It was here decided to fit the critical temperature and the critical pressure of the pseudocomponent C 11 -C 19 on all the available experimental data except the MMP and then to predict the MMP with the fitted parameters. Such a procedure is here justified because the partial volumes during MCT are poorly predicted and because the saturation pressures during the swelling tests are not so good. We hope that by reducing the deviations on these data, we will improve the prediction of the MMP. It is here recalled that the question we (d-f) Swelling test performed on fluid F5 (gas G5b is injected). The average deviations of the saturation pressures, swollen volumes, and densities of the swollen oil are respectively equal to 4.75% (d), 0.69% (e), and 1.96% (f).

Table 2. Molar Composition of the 13 Crude Oils and of the 17 Injected Gases Used in This Paper

want to answer is, is it still necessary to measure the MMP? In other words, may other gas injection experiments such as swelling test or MCT substitute for slim tube test. The results obtained for the first three systems (F5 + G5b, F6 + G6b, F12 + G12) are given in Table 4 and are illustrated in Figure 2. From this table it becomes obvious that for these three examples the MMP may be estimated with high accuracy (average deviation of 3.8%) after fitting only two parameters on data including swelling test and MCT. The best result (1.1% deviation) is obtained for the system F5 + G5b in which four contacts were performed during the MCT. Moreover the worst result (5.6% deviation) is obtained for system F6 + G6b in which no MCT experiments are available. An intermediate result (4.7% deviation) is obtained for system F12 + G12 in which two experiments of backward MCT were conducted but for which only one contact was performed. These results seem to indicate that to predict accurately the MMP, it is necessary to perform a swelling test but above all MCT experiments. The accuracy on the MMP seems to be directly proportional to the contact number (four con-tacts are however enough to obtain the MMP with a deviation of 1%).

Concerning the last system F13 + G13, it was not possible using the same strategy to obtain accurate results. This is because the injected gas is pure CO 2 .It is indeed well-known that the PR EOS is not able to properly predict the phase behavior of mixtures containing a large amount of CO 2 . As an example, using the compositional model described in this paper, it is impossible to compute a dew point pressure as soon as the amount of CO 2 injected in the reservoir fluid is higher than 85%. Indeed, the saturation pressure becomes higher than 10 6 bar, which is the maximum working pressure of our software. This phenomenon is visible in Figure 3. A similar behavior was observed during the simulation of a swelling test when pure nitrogen (N 2 ) or pure hydrogen sulfide (H 2 S) is injected. In this last particular case, it was thus decided to include the MMP in the fitting procedure. The results are given in Table 5. This table provides evidence that by fitting the critical pressure and the critical temperature of the C 11 -C 19 pseudocomponent, it was possible to exactly match the experimental MMP value. However some other data, such as the bubble point reservoir pressure, are not reproduced very accurately. What is also here very interesting to outline is that the partial volumes during the MCT experiments are also very accurately calculated with the fitted parameters. Once more, there is a correlation between a good prediction of these relative volumes and the MMP. However because of the behavior shown in Figure 3, it was not possible to predict accurately the MMP after fitting two parameters on data which did not include the MMP.

Conclusion

In the first part of this paper, a slightly modified version of a previously developed compositional model coupled to a predictive cubic EOS is tested. Such a model is able to predict the classical PVT data (relative volumes, bubble pressures, swelling test data, . . .) with high accuracy. However, the partial volumes during a MCT and the MMP are not always precisely predicted. The classical solution for this kind of problem is to fit some parameters of the EOS on the experimental data in order to reduce the deviations between experimental and calculated values. When no other solutions were available, such a procedure was used with high success in this paper. In particular our software is uniquely able to take into account the MMP described by a mixed condensing/vaporizing mechanism in the fitting procedure. The main problem with such a procedure comes from the MMP needs to be experimentally determined. But because the determination of the MMP with a slim tube is a very expensive experiment, the petroleum companies wish to avoid such an experiment. This paper gives a solution to this problem. Indeed, we have shown that when the injected gas is not pure CO 2 (and probably not pure N 2 or pure H 2 S), it is enough to fit only two parameters of the EOS on data including classical PVT data + swelling data + MCT data and then to predict the MMP. During MCT experiments, this paper shows that it is very important to perform a few contacts (four contacts seems to be enough). We can then conclude that the slim tube test may be replaced by swelling tests and MCTs, which are much cheaper.

In conclusion, when the injected gas is not pure CO 2 (or N 2 or H 2 S), it is not necessary to measure the MMP.

Figure 1 .

 1 Figure 1. Illustration of the predictive power of the proposed compositional model. Symbols + and 0 are the experimental data points. The continuous lines are the predicted values. (a) Constant mass expansion of fluid F5 at T/K ) 394.25. The average deviation between calculated and experimental relative volumes is equal to 1.03%. (b) Differential vaporization of fluid F11 at T/K ) 373.75. The average deviation between calculated and experimental relative volumes is equal to 1.57%. (c) (P,T) phase envelope of fluid F12 for which four bubble point pressures were measured. The average deviation of these four pressures is equal to 2.92%. (d-f) Swelling test performed on fluid F5 (gas G5b is injected). The average deviations of the saturation pressures, swollen volumes, and densities of the swollen oil are respectively equal to 4.75% (d), 0.69% (e), and 1.96% (f).

Figure 2 .

 2 Figure 2. Illustration of the results obtained when two parameters are fitted on all the available experimental data points except the MMP. Symbols + and 0 are the experimental data points. The continuous lines are the calculated values with the fitted parameters. (a) Constant mass expansion of fluid F6 at T/K ) 383.15. The average deviation between calculated and experimental relative volumes is equal to 0.71%. (b-d) Swelling test performed on fluid F6 (gas G6b is injected). The average deviations of the saturation pressures, swollen volumes, and densities of the swollen oil are respectively equal to 1.91% (b), 0.42% (c), and 3.41% (d).

Figure 3 .

 3 Figure 3. Swelling test performed on fluid F13 (pure CO2 is injected): evolution of the saturation pressure versus quantity of injected gas. Experimental data points are shown as boxes. The continuous lines are the predicted values (the calculation becomes impossible as soon as xCO 2 becomes higher than 0.85).

Table 1 . Results Obtained with the Purely Predictive Compositional Model Used in This Study

 1

						crude	stock
						oil	tank					isothermal and
		constant mass	differential	bubble	oil					isobaric multicontact	minimum miscibility
		expansion	vaporization	point	density	isothermal swelling test		test			pressure
	fluid		∆V rel		∆Vrel	∆Psat	∆F		∆Psat.	∆Vswell	∆Fsat.			∆Pvol		mech-	∆MMP
	name T/K	(%)	T/K	(%)	(%)	(%)	gas	(%)	(%)	(%)	gas	type	(%)	gas	anism	(%)
	F1	374.85 1.43 (9) a			-5.45	1.79								G1	C/V	-2.4
	F2	372.05 1.54 (11)			7.55									G2	C/V	2.8
	F3	387.35 0.35 (9) 387.35 2.07 (11) a -1.12	0.44	G3 (6) a 2.16	1.26	2.42				G3	C/V	-2.5
	F4	388.15 0.91 (11)			3.61	3.40	G4 (7)	3.03	0.26	0.34				G4	b	2.3
	F5	394.25 1.03 (18) 394.25 0.51 (12)	4.39	0.36	G5a (6) 5.05	0.37	1.01				G5a C/V	-2.8
								G5b (7) 4.75	0.69	1.96 G5c	c	37.7 G5d C/V	9.4
	F6	383.15 1.15 (12)			-4.58	1.60								G6a C/V	-3.5
								G6b (6) 9.05	0.25	1.16				G6b C/V	-9.8
	F7	393.15 1.38 (12)			6.02	4.01								G7	C/V	-5.4
	F8	393.15 0.49 (1)			-1.07	1.25								G8	C/V	-21.3
	F9	394.25 0.20 (10) 394.25 2.37 (12)	0.57	2.21								G9	b	-5.4
		303.25 0.67 (12)			6.12
		352.45 0.39 (10)			2.19
		417.85 0.37 (13)			-0.49
	F10 394.25 2.01 (9)			-7.14	0.52								G10 C/V	14.1
	F11 373.75 0.72 (9) 373.75 1.57 (11) -2.88	1.29								G11 C/V	-16.0
	F12 376.45 0.49 (9) 376.45 2.13 (11)	3.14	0.35	G12 (4) 4.49	0.25	2.86 G12	d	-3.6 G12 C/V	8.6
		299.85 0.44 (9)			5.22
		360.95 0.31 (9)			1.70							d	43.9
		394.25 0.41 (9)			1.66
	F13 377.55 0.80 (15) 377.55 0.72 (14) -5.09	0.52	G13 (6) 9.54	1.97	1.20 G13	c	5.0 G13 C/V	7.4
													e	32.5
	av		0.83		1.52	3.68	1.48		5.41	0.73	1.48			24.3			7.6
			(197) a		(71)	(19)	(12)		(42)	(42)	(42)			(13)			(15)

a

The number in parentheses is the number of data points. b Vaporizing. c Reverse (four contacts). d Reverse (one contact). e Forward (three contacts).

Table 3 . Deviations Obtained after Fitting Only Two Parameters (TC,C 11-C19 and PC,C 11-C19 ) on All the Available Experimental Data (the MMP is Included in the Fitting Procedure)

 3 The number in parentheses is the average deviation obtained with the predictive model.

		constant mass	differential	crude oil	stock tank	minimum miscibility
		expansion,	vaporization,	bubble point,	oil density,	pressure,
	fluid name	∆Vrel (%)	∆Vrel (%)	∆Psat. (%)	∆F (%)	∆MMP (%)
	F7	1.09		3.64	-3.40	-0.4
	F8	0.38		4.84	-2.55	-0.5
	F9	0.41	2.86	1.44	-2.21	-0.6
		0.67		6.12
		0.45		2.57
		0.46		0.72
	F10	1.55	-	-5.88	8.10	7.5
	F11	0.58	1.78	2.73	0.09	-0.8
	av	0.73 (0.79) a	2.34 (1.99)	3.49 (3.31)	3.27 (1.86)	2.0 (12.4)

a

Table 4 . Deviations Obtained after Fitting Only Two Parameters on All the Available Experimental Data Except the MMP and MMP Prediction with the Fitted Parameters

 4

				crude	stock
		constant		oil	tank
		mass	differential	bubble	oil					isothermal and	minimum miscibility
		expansion	vaporization	point	density	isothermal swelling test	isobaric MCT		pressure a
	fluid		∆Vrel	∆Psat.			∆Psat	∆Vswell	∆Fsat.			∆Pvol		mech-	∆MMP
	name	∆Vrel (%)	(%)	(%)	∆F (%)	gas	(%)	(%)	(%)	gas type	(%)	gas	anism	(%)
	F5	0.88	0.44	0.56	1.14	G5b	1.26	0.27	1.54	G5b	b	16.5	G5b	C/V	1.1
	F6	0.71		3.50	2.52	G6b	1.91	0.42	3.41				G6b	C/V	5.6
	F12	0.33	1.44	-0.90	3.43	G12	1.70	0.21	3.45	G12	c	-6.7	G12	C/V	-4.7
		0.09		1.84
		0.69		-2.24							c	15.5
		0.74		-2.67
	av	0.62	0.92	1.95	2.36		1.59	0.31	2.65			14.7			3.8
		(0.72) d	(1.28)	(3.45)	(0.77)		(6.21) (0.43) (1.89)			(33.1)			(9.3)

a Prediction. b Reverse (four contacts). c Reverse (one contact). d The number in parentheses is the average obtained with the predictive model.

Table 5 . Deviations Obtained after Fitting Only Two Parameters on All the Available Experimental Data Including the MMP

 5

				crude	stock
		constant		oil	tank
		mass	differential	bubble	oil					isothermal and	minimum miscibility
		expansion	vaporization	point	density	isothermal swelling test	isobaric MCT		pressure
	fluid	∆Vrel	∆Vrel	∆Psat.	∆F		∆Psat.	∆Vswell	∆Fsat			∆Pvol		mech-	∆MMP
	name	(%)	(%)	(%)	(%)	gas	(%)	(%)	(%)	gas type	(%)	gas	anism	(%)
	F13	1.82	1.37	-10.1	-8.10	G13	6.17	2.66	5.56 G13	a	4.5	G13	C/V	0.3
											b	2.7

a Reverse (four contacts). b Forward (three contacts).

Fax number: (+33) 3 83.17.51.52. Phone number: (+33) 3 83.17.50.81. E-mail: Jean-Noel.Jaubert@ensic.inpl-nancy.fr.† Institut National Polytechnique de Lorraine. ‡ Parc Technologique de Pre ´Milliet.

Acknowledgment

This work was sponsored by the French petroleum company TOTALFINAELF. The authors wish to especially thank Dr. Marc Durandeau (expert in thermodynamics) and Cyril Le Goff (head of PVT section) for helpful discussions.


Sibo Cheng
email: sibo.cheng@imperial.ac.uk

César Quilodrán-Casas

Rossella Arcucci
email: r.arcucci@imperial.ac.uk

Reduced order surrogate modelling and Latent Assimilation for dynamical systems

Keywords: Reduced order modelling, Deep learning, Data assimilation, Adversarial neural networks

For high-dimensional dynamical systems, running high-fidelity physical simulations can be computationally expensive. Much research effort has been devoted to develop efficient algorithms which can predict the dynamics in a low-dimensional reduced space. In this paper, we developed a modular approach which makes use of different reducedorder modelling for data compression. Machine learning methods are then carried out in the reduced space to learn the dynamics of the physical systems. Furthermore, with the help of data assimilation, the proposed modular approach can also incorporate observations to perform real-time corrections with a low computational cost. In the present work, we applied this modular approach to the forecasting of wildfire, air pollution and fluid dynamics. Using the machine learning surrogate model instead of physics-based simulations will speed up the forecast process towards a real-time solution while keeping the prediction accuracy. The data-driven algorithm schemes introduced in this work can be easily applied/extended to other dynamical systems.

Introduction

Reduced-order modelling is a reduced dimensionality model surrogate of an existing system in a reduced-order space. reduced-order modelling (ROM), in combination with machine learning (ML) algorithms is of increasing interest of research in engineering and environmental science. This approach improves the computational efficiency for high-dimensional systems. Since forecasting the full physical space is computationally costly, much effort has been given to develop ML-based surrogate models in the pre-trained reduced-order space.

In recent years, the algorithm schemes which combine ROM and ML surrogate models have been applied in a variety of engineering problems, including computational fluid dynamics [START_REF] Casas | A Reduced Order Deep Data Assimilation model[END_REF][START_REF] Quilodrán-Casas | Adversarial autoencoders and adversarial LSTM for improved forecasts of urban air pollution simulations[END_REF], numerical weather prediction [START_REF] Knol | Deep learning for solar irradiance nowcasting: A comparison of a recurrent neural network and two traditional methods[END_REF] and nuclear science [START_REF] Phillips | An autoencoder-based reduced-order model for eigenvalue problems with application to neutron diffusion[END_REF], among others, to speed up computational models without losing the resolution and accuracy of the original model. Typically, the first stage consists of reducing the dimension of the problem by compression methods such as Principal Component Analysis (PCA), autoencoder (AE), or a combination of both [START_REF] Casas | Urban air pollution forecasts generated from latent space representations[END_REF][START_REF] Phillips | An autoencoder-based reduced-order model for eigenvalue problems with application to neutron diffusion[END_REF]. Solutions from the original computational models (known as snapshots) are then projected onto the lower-dimensional space, and the resulting snapshot coefficients are interpolated in some way, to approximate the evolution of the model.

In order to incorporate real-times observations, data assimilation (DA), originally developed in meteorological science, is a reference method for system updating and monitoring. Some recent researches [START_REF] Cheng | Observation data compression for variational assimilation of dynamical systems[END_REF][START_REF] Amendola | Data assimilation in the latent space of a neural network[END_REF][START_REF] Casas | A reduced order deep data assimilation model[END_REF][START_REF] Liu | Enkf data-driven reduced order assimilation system[END_REF] have focused on combining DA algorithms and ROMs so that the system correction/adjusting can be performed with a low computational cost. Adversarial training and Generative Adversarial Networks (GAN), introduced by [START_REF] Goodfellow | Generative adversarial nets[END_REF], has also been used with ROM. Data-driven modelling of nonlinear fluid flows incorporating adversarial networks has been successfully studied previously [START_REF] Cheng | Data-driven modelling of nonlinear spatio-temporal fluid flows using a deep convolutional generative adversarial network[END_REF]. GANs are also being used to capture physics of molecular dynamics [START_REF] Wu | Deep generative markov state models[END_REF] and have potential to aid in the modeling and simulation of turbulence [START_REF] Kim | Deep fluids: A generative network for parameterized fluid simulations[END_REF].

The aim of this work is to create general workflows in order to tackle different applications combining DA and ML approaches. In the present work, we propose a modular approach, which combines ROM, ML surrogate models, and DA for complex dynamic systems with applications in computational fluid dynamic (CFD), wildfire spread and air pollution forecasting. The algorithms described in this work can be easily applied/extended to other dynamical systems. Numerical results in both applications show that the proposed approach is capable of real-time predictions, yielding an accurate result and considerable speed-up when compared to the computational time of simulations.

The paper is structured as follows: Section 2 presents the modular approach and its components. The applications are shown in Section 3 and 4. And finally, Discussions and Conclusions are presented in Section 5.

Components of the modular approach

The modular approach shown in this paper can be summarised by Figure 1. The state model (u t ) is compressed using ROMs approaches such as PCA, AE or a combination of both, followed by a ML-based forecast in the reduced space. This forecast is then corrected using DA, incorporating real-time observations (v t ). This is an iterative process that can be used to improve the starting point of the next time-level forecast, thus improving its accuracy [START_REF] Asch | Data assimilation: methods, algorithms, and applications[END_REF].

Reduced order modelling

In this section, we introduce two types of ROMs, namely the PCA and the convolutional autoencoder (CAE). Fig. 1: Flowchart of the modular approach, modified using figure 1 of [START_REF] Buizza | Data learning: Integrating data assimilation and machine learning[END_REF] 2.1.1 Principle component analysis principal component analysis, also known as Karhunen-Loève transform or Hotelling transform, is a reference method of ROM via an orthogonal linear projection. This approach has been widely applied in dynamical systems [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. II. Symmetries and transformations[END_REF] with snapshots at different time steps. Applications can be found in a large range of engineering problems, including numerical weather prediction [START_REF] Jaruszewicz | Application of pca method to weather prediction task[END_REF], hydrology [START_REF] Cheng | Observation data compression for variational assimilation of dynamical systems[END_REF] or nuclear engineering [START_REF] Gong | Data-enabled physics-informed machine learning for reduced-order modeling digital twin: Application to nuclear reactor physics[END_REF]. More precisely, a set of n u simulated or observed fields {u t0,t1,..tn u -1 } at different time are flattened and combined vertically to a matrix,

U = u t0 u t1 ... u tn u -1 . (1)
The principle components are then extracted by computing the empirical covariance matrix, that is,

C u = 1 n u -1 UU T = L U D U L U T , (2)
where each column of L U which represents an eigenvector of C u and D U is the associated eigenvalue diagonal matrix. The dynamical field u t can then be compressed to

ũt = L U,q T u t , (3)
where ũt denotes the compressed state vectors and q is the truncation parameter and the reconstruction to the full physical space reads

u PCA t = L U,q ũt = L U,q L U,q T u t . (4)
2.1.2 Convolutional autoencoder PCA, by design, is a linear method for ROM. In the last several years, much effort has been given to dimension reduction for chaotic dynamical systems via deep learning (DL) AEs [START_REF] Dong | A review of the autoencoder and its variants: A comparative perspective from target recognition in synthetic-aperture radar images[END_REF]. Typically, an AE is an unsupervised neural network (NN) which consists two parts: an encoder E which maps the input variables to latent (i.e., compressed) vectors and a decoder D which performs reconstructions from the low-dimensional latent space to the full physical space. These processes can be summarised as:

ũt = E(u t ) and u AE t = D(ũ t ). (5)
It is found that employing convolutional layers in AEs is helpful to i) reduce the number of parameters for high-dimensional systems, and ii) take into account local spatial patterns for structured data (e.g., images and times series) [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF]. Following this idea, CAE was developed [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF][START_REF] Rawat | Deep convolutional neural networks for image classification: A comprehensive review[END_REF] where both the encoder E and the decoder D consist of a series of convolutional layers.

In general, the encoder and the decoder of an AE should be trained jointly, for instance, with the Mean square error (MSE) or the Mean absolute error (MAE) loss function of reconstruction accuracy, i.e.,

J MSE (E, D) = nu-1 j=0 ||u tj -D • E(u tj )|| 2 n u , J MAE (E, D) = nu-1 j=0 ||u tj -D • E(u tj )|| n u . (6)
2.2 Machine Learning for surrogate models

Recurrent neural network

Long short-term memory networks neural networks, introduced in [START_REF] Hochreiter | Long short-term memory[END_REF], is a variant of recurrent neural network (RNN), capable of dealing long term dependency, and vanishing gradient problems that traditional RNN could not handle. One of the components of our modular approach in the present work, is the sequenceto-sequence long short-term memory networks (LSTM). The LSTM learns the dynamics in the latent space from compressed training data. LSTM can be unidirectional or bidrectional. Recently developed Bidirectional LSTM (BDLSTM) [START_REF] Schuster | Bidirectional recurrent neural networks[END_REF] differs from the unidirectional ones, as the latter can capture the forward and backward temporal dependencies in spatiotemporal data [START_REF] Cui | Stacked bidirectional and unidirectional LSTM recurrent neural network for forecasting network-wide traffic state with missing values[END_REF][START_REF] Huang | Bidirectional lstm-crf models for sequence tagging[END_REF][START_REF] Liu | Bidirectional lstm with attention mechanism and convolutional layer for text classification[END_REF][START_REF] Quilodrán-Casas | Digital twins based on bidirectional lstm and gan for modelling the covid-19 pandemic[END_REF] . LSTMs are widely recognised as one of the most effective sequential models for times series predictions in engineering problems [START_REF] Tekin | Spatiotemporal weather forecasting and attention mechanism on convolutional lstms[END_REF][START_REF] Nakamura | Convolutional neural network and long short-term memory based reduced order surrogate for minimal turbulent channel flow[END_REF].

The LSTM network comprises three gates: input (i t k ), forget (f t k ), and output (o t k ); a block input, a single cell c t k , and an output activation function. This network is recurrently connected back to the input and the three gates. Due to the gated structured and the forget state, the LSTM is an effective and scalable model that can deal with long-term dependencies [START_REF] Hochreiter | Long short-term memory[END_REF]. The vector equations for a LSTM layer are:

i t k = ϕ(W xi u t k + W Hi H t k-1 + b i ) f t k = ϕ(W xf u t k + W Hf H t k-1 + b f ) o t k = ϕ(W xo u t k + W Ho H t k-1 + b o ) c t k = f t k • c t k-1 + i t k • tanh(W xc u t k + W Hc H t k-1 + b c ) H t k = o t k • tanh(c t k ) (7)
where ϕ is the sigmoid function, W are the weights, b i,f,o,c are the biases for the input, forget, output gate and the cell, respectively, x t k is the layer input, H t k is the layer output and • denotes the entry-wise multiplication of two vectors. This is the output of a unidirectional LSTM.

For a BDLSTM, the output layer generates an output vector u t k :

u t k = ψ( --→ H t k , ← -- H t k ) ( 8
)
where ψ is a concatenating function that combines the two output sequences, forwards and backwards, denoted by a right and left arrow, respectively.

Adversarial network

The work of [START_REF] Goodfellow | Generative adversarial nets[END_REF] introduced the idea of adversarial training and adversarial losses which can also be applied to supervised scenarios and have advanced the state-of-the-art in many fields over the past years. Additionally, robustness may be achieved by detecting and rejecting adversarial examples by using adversarial training [START_REF] Shafahi | Adversarial training for free! In[END_REF]. GAN are a network trained adversarially. The basic idea of GAN is to simultaneously train a discriminator and a generator, where the discriminator aims to distinguish between real samples and generated samples. By learning and matching the distribution that fits the training data x, the aim is that new samples, sampled from the matched distribution formed by the generator, will produce, or generate, 'realistic' features from the latent vector z The GAN is composed by a discriminator network (D) and a generator network (G) The GAN losses, binary cross-entropy, therefore, can be written as:

L adv D = - logD(x)) + log(1 -D(G(z)) L adv G = - log(D(G(z)) (9)
This idea can be developed further if we consider similar elements of the adversarial training of GAN and applied to other domains, e.g. time-series, extreme events detection, adversarial attacks, among others.

Data assimilation

Data assimilation algorithms aim to estimate the state variable u relying on a prior approximation u b (also known as the background state) and a vector of observed states v. The theoretical value of the state vector is denoted by u true , so called the true state, which is out of reach in real engineering problems. Both the background and the observation vectors are supposed to be noisy in DA, characterised by the associated error covariance matrices B and R, respectively, i.e.,

B = Cov(ϵ b , ϵ b ), R = Cov(ϵ o , ϵ o ), (10)
with the prior errors ϵ b and ϵ o defined as:

ϵ b = u b -u true ϵ o = H(u true ) -v. (11)
Since the the true states are out of reach in real applications, the covariance matrices B and R are often approximated though statistical estimations [START_REF] Cheng | Background error covariance iterative updating with invariant observation measures for data assimilation[END_REF][START_REF] Desroziers | Diagnosis and adaptive tuning of observation-error parameters in a variational assimilation[END_REF]. The H function in equation 11 is called the transformation operator, which maps the state variables to the observable quantities. H(u true ) is also known as the model equivalent of observations.

By minimizing a cost function J defined as

J(u) = 1 2 (u -u b ) T B -1 (u -u b ) + 1 2 (v -H(v)) T R -1 (v -H(v)) = J b (u) + J o (u) (12)
DA approaches attempt to find an optimally weighted a priori analysis state,

u a = argmin u J(u) . (13)
The B and the R matrices, determining the weights of background and observation information (as shown in equation 12), is crucial in DA algorithms [START_REF] Tandeo | A review of innovation-based methods to jointly estimate model and observation error covariance matrices in ensemble data assimilation[END_REF][START_REF] Cheng | Observation error covariance specification in dynamical systems for data assimilation using recurrent neural networks[END_REF]. When H can be approximated by a linear function H and the error covariances B and R are well specified, equation 12 can be solved via Best Linear Unbiased Estimator (BLUE) [START_REF] Cheng | Background error covariance iterative updating with invariant observation measures for data assimilation[END_REF]:

u a = u b + K(v -Hu b ) (14)
where K denotes the Kalman gain matrix,

K = BH T (HBH T + R) -1 . (15)
The optimisation of equation 12 often involves gradient descent algorithms (such as "L-BFGS-B") and adjoint-based numerical techniques. In the proposed modular approach of the present paper, we aim to perform DA in the low-dimensional latent space to reduce the computational cost, enabling a real-time model updating. The latent assimilation (LA) approach was first introduced in the work of [START_REF] Amendola | Data assimilation in the latent space of a neural network[END_REF] for CO 2 spread modeling. A generalised Latent Assimilation algorithm was proposed in the recent work of [START_REF] Cheng | Generalised latent assimilation in heterogeneous reduced spaces with machine learning surrogate models[END_REF]. The observation quantities v t are first preprocessed to fit the space of the state variables u t , i.e.,

ṽt = E(v t ). (16)
As a consequence, the transformation operator becomes the identity function in the latent space, leading to the loss function of LA:

J( ũt ) = 1 2 (ũ t -ũt,b ) T B -1 (ũ t -ũt,b ) + 1 2 (ũ t -ṽt ) T R -1 (ũ t -ṽt ), (17)
where the latent background state ũt,b is issued from the RNN predictions as mentioned in section 2.2.1. The analysis state,

u t,a = argmin ut J(ũ t ) , (18)
can then replace the background prediction ũt,b , which can be used as the starting-point for the next-level prediction in ML algorithms.

Application to Wildfires

The first application is real-time forecasting of wildfire dynamics. Wildfires have increasing attention recently in fire safety science world-widely, and it is an extremely challenging task due to the complexities of the physical models and the geographical features. Real-time forecasting of wildfire dynamics which raises increasing attention recently in fire safety science world-widely, is extremely challenging due to the complexities of the physical models and the number of geographical features. Running physics-based simulations for large-scale wildfires can be computationally difficult, if not infeasible. We applied the proposed modular approach for fire forecasting in near real-time, which combines reducedorder modelling, recurrent neural networks RNN, data assimilation DA and error covariance tuning. More precisely, based on snapshots of dynamical fire simulations, we first construct a low-dimensional latent space via proper orthogonal decomposition or convolution AE. A LSTM is then used to build sequence-tosequence predictions following the simulation results projected/encoded in the reduced space. In order to adjust the prediction of burned areas, latent DA coupled with an error covariance tuning algorithm is performed with the help of daily observed satellite wildfire images as observation data. The proposed method was tested on two recent large fire events in California, namely the Buck fire and the Pier fire, both taking place in 2017 as illustrated in figure 2. We first employed an operational cellular automata (CA) fire spread model [START_REF] Alexandridis | A cellular automata model for forest fire spread prediction: The case of the wildfire that swept through spetses island in 1990[END_REF] to generate the training dataset for ROM and RNN surrogate modelling. This CA model is a probabilistic simulator which takes into account a number of local geophysical features, such as vegetation density (see figure 2) and ground elevation. Once the latent space is acquired, the ML-based surrogate model is then trained using the results of stochastic CA simulations in the corresponding area of fire events. With a much shorter online execution time, the so-obtained data-driven model provides similar results as physics-based CA simulations (stochastic) in the sense that the mean and the standard deviation of CA-CA and CA-LSTM differences are similar as shown in figure 3 for the Pier fire. In fact, the ROMand ML-based approach run roughly 1000 times faster than the original CA model as shown in figure 3(b). Each step in CA-LSTM predictions is roughly equivalent to 30 minutes in real time while the satellite observations are of daily basis. The latter is used to adjust the fire prediction consistently since the actual fire spread also depends heavily on other factors such as real-time climate or human interactions which are not included in the CA modelling. The evolution of the averaged relative root mean square error (R-RMSE) is shown in figure 4 . The numerical results show that, with the help of error covariance tuning [START_REF] Desroziers | Diagnosis and adaptive tuning of observation-error parameters in a variational assimilation[END_REF][START_REF] Cheng | Error covariance tuning in variational data assimilation: application to an operating hydrological model[END_REF], DA manages to improve the model prediction accuracy in both fire events.

Application to Computational Fluid Dynamics and air pollution

Similar to the wildfire problem, we also present an general workflow to generate and improve the forecast of model surrogates of CFD simulations using deep learning, and most specifically adversarial training. This adversarial approach aims to reduce the divergence of the forecasts from the underlying physical model. Our two-step method, similar to the wildfire application, integrates a PCA AE with adversarial LSTM networks. Once the reduced-order model (ROM) of the CFD solution is obtained via PCA, an adversarial autoencoder (AAE) is used on the principal components time series. Subsequently, a LSTM model is adversarially trained, named adversarial LSTM (ALSTM), on the latent space produced by the principal component adversarial autoencoder (PC-AAE) to make forecasts. Here we show, that the application of adversarial training improves the rollout of the latent space predictions. UK. The 3D case (720m × 676m × 250m) is composed of an unstructured mesh including m = 100, 040 nodes per dimension and n = 1, 000 time-steps.

The two-dimensional (2D) CFD case study was performed using Thetis [START_REF] Kärnä | Thetis coastal ocean model: discontinuous galerkin discretization for the threedimensional hydrostatic equations[END_REF] and the three-dimensional (3D) CFD simulations were carried out using Fluidity [START_REF] Davies | Fluidity: A fully unstructured anisotropic adaptive mesh computational modeling framework for geodynamics[END_REF]. For these domains, the framework was trained and validated on the first 1000 time-steps of the simulation, and tested on the following 500 time-steps.

A PCA as applied to a 2-dimensional velocity field (ms -1 ) in Flow past the cylinder and likewise to the velocities of the 3D model. The full-rank PCs were used as input for the AAE and divided in 3 different experiments named LS τ and compared to the corresponding reconstruction x τ with τ = {2, 4, 8} PCs. The results of the mean absolute error using the different dimension reduction approaches are shown in Figure 6a for the flow past the cylinder case. The AAE outperforms a simple truncation of the PCs in both domains.

In terms of forecasting, our framework generalises well on unseen data 6b. This is because of the Gaussian latent space obtained with the adversarial AE constraints the further predictions and forces the predictions back into the distribution. Furthermore, the adversarial training of the LSTM learns how to stay within the distribution data. Followed by the training of an adversarial LSTM, we can assess the forecasts using our workflow. An ensemble of 50 different starting points from the test dataset were used to be forecasted in time for 100 timelevels. The ensemble of Mean Absolute Errors results are based on a dimension reduction 8 dimensions in the latent space of the AAE, which is a compression of 5 orders of magnitude. The error percentage of the means of these forecasts is 5% in the test dataset.

Conclusions

In the present paper, we introduced a ROM-and ML-based modular approach for efficient predictions of high-dimensional dynamical systems. In addition, this method can also incorporate real-time observations for model correction/adjusting with a low computational cost. A variety of ROM and RNN approaches can be included in the algorithm scheme regarding different applications. The replacement of the physics-based simulation/resolution by these models will speed up the forecast process towards a real-time solution. And, the application of adversarial training could potentially produce more physically realistic scenarios. We showed the strength of the proposed method in predicting wildfire spread and air pollution diffusion in this work. Furthermore, this framework is data-agnostic and could be applied to different physical models when enough data is available.

Fig. 2 :

 2 Fig. 2: Observed burned area of the Buck fire and the Pier fire at the first day

Fig. 3 :

 3 Fig. 3: Model difference and computational time

Fig. 4 :Fig. 5 :

 45 Fig. 4: Relative error of pure LSTM and assimilated prediction compared to satellite observations

Fig. 6 :

 6 Fig. 6: Top: Mean absolute error (MAE) of the different compression methods on the velocities (ms -1 ) of the flow past the cylinder domain. Red, magenta, and cyan represents a compression to 8, 4, and 2 dimensions, respectively. Circle and triangle markers are a PC and a PC-AAE compression, respectively. Bottom: Ensemble, of 50 different starting points, of the MAE of the forecast, for 100 time-levels, of the velocities (ms -1 ) of the 3D Urban Air Pollution domain. Blue is the original test data, and orange is the prediction obtained with the PC-AAE and the ALSTM.

Acknowledgements

This research is funded by the Leverhulme Centre for Wildfires, Environment and Society through the Leverhulme Trust, grant number RC-2018-023. This work is partially supported by the EP/T000414/1 PREdictive Modelling with

QuantIfication of UncERtainty for MultiphasE Systems (PREMIERE) and by the EPSRC grant EP/T003189/1 Health assessment across biological length scales for personal pollution exposure and its mitigation (INHALE).