Yingpu Xie

Kuo Zeng
email: zengkuo666@hust.edu.cn

Gilles Flamant

Haiping Yang

Nian Liu

Xiao He

Xinyi Yang

Ange Nzihou

Hanping Chen

Solar pyrolysis of cotton stalk in molten salt for bio-fuel production

Keywords:

Solar pyrolysis of cotton stalk was carried out in a molten salt reactor heated by 4 kW solar simulator. The effects of pyrolysis temperature and mass ratio of molten salt to biomass on pyrolysis products properties were investigated. The use of molten salt as pyrolysis media increased gas yield. At 850 C, the gas yield (mainly of CO and H 2 ) continued to rise from 41.35 wt% to 82.57 wt% when mass ratio of molten salt to cotton stalk increased from 0 to 10. Pyrolysis in molten salt significantly decreased bio-oil acids and phenols, while increased aromatics among pyrolysis temperature range of 450e850 C. There was positive correlation between the increased content of aromatics and mass ratio of molten salt to cotton stalk (from 0.5 to 10). The bio-char carbon content showed a general decreasing trend while oxygen, BET surface area and pore volume increased with using molten salt as pyrolysis media. Bio-char obtained from CS1MS5 pyrolysis at 850 C had the highest BET surface area of 972.57 m 2 /g and the biggest total pore volume of 0.6203 cm 3 /g. High quality pyrolysis products with more uniform chemistry suggest catalytic reactions occur inside the solar reactor due to the intermediates degradation with molten salt.

Introduction

Due to the negative environmental impacts of fossil fuels and the increasing global energy demands, biofuels are receiving increasing attention as the best short-term substitute for petroleum fuels. Pyrolysis is a promising thermochemical process for converting biomass into biofuels (bio-gas, bio-oil and bio-char) [START_REF] Neves | Characterization and prediction of biomass pyrolysis products[END_REF]. However, the thermal energy driving the pyrolysis process is typically drawn from combustion of the feedstock or other fossil energy sources. To produce a truly carbon-neutral fuel and maximize the yield of biofuel per unit biomass consumed, the energy driving the pyrolysis reactions must come from a renewable energy source. Concentrated solar energy fits the requirements of a heat source for an externally heated pyrolysis process. Solar radiation is concentrated using mirrors to provide temperatures high enough to drive biomass pyrolysis reactions that are able to store this solar thermal energy in a chemical form as biofuels [2e4].

Supplying process heat with concentrated solar energy has numerous advantages over conventional pyrolysis [START_REF] Zeng | Solar pyrolysis of carbonaceous feedstocks: a review[END_REF]. There is no need for partial combustion which allowing all of the feedstock to be upgraded into bio-fuel; the product gas contains low CO 2 and N 2 , and no other combustion by-products; the energetic value of the product is greater than that of the feedstock. Besides, pyrolysis driven by direct concentrated solar irradiation has the benefits of high temperatures and rapid heating rates [START_REF] Lede | Solar thermochemical conversion of biomass[END_REF]. The biomass can be rapidly heated to around 500 C and thus maximizes exposure to the higher temperatures that favor vapors formation. Hence, high yields of bio-oil will be obtained at normal pyrolysis temperatures of about 500 C because vapors passing from the high temperature region (the focus of concentrated solar radiation) to the low temperature zone (the surrounding gas) are quickly condensed, thus minimizing the secondary decomposition of primary tar into gas [START_REF] Hopkins | The flash pyrolysis of cellulosic materials using concentrated visible light[END_REF]. Fuel calorific value of carbonaceous feedstock can be upgraded by direct solar pyrolysis process, which enables the storage of solar energy in chemical form as bio-oil, bio-char and bio-gas whose proportion mainly depending on the temperature [START_REF] Zeng | Modeling of beech wood pellet pyrolysis under concentrated solar radiation[END_REF][START_REF] Zeng | Combined effects of initial water content and heating parameters on solar pyrolysis of beech wood[END_REF]. The previous work on pyrolysis driven by direct irradiation of concentrated solar energy has demonstrated many benefits, however issues related to transient sensitivity, poor heat transfer and tar production have been persistent motivating the search for an improved approach to solar pyrolysis. Concentrated solar energy can be effectively stored in molten salt, then providing high temperature heat for biomass pyrolysis [START_REF] Nygård | Review of thermal processing of biomass and waste in molten salts for production of renewable fuels and chemicals[END_REF]. The use of molten salts allowed the continuous operation of solar pyrolysis system. Besides, the biomass particles can be dispersed in molten salt, which is beneficial to enhancing the bio-fuels production rate and yield. This can be attributed to two main reasons [START_REF] Hathaway | Solar gasification of biomass: kinetics of pyrolysis and steam gasification in molten salt[END_REF]: fast heat transfer rate and catalytic effect of alkali metal-based catalyst. The fast heat transfer was attributed to a high convective heat transfer with high heat capacity [START_REF] Hathaway | Steam gasification of plant biomass using molten carbonate salts[END_REF]. The alkali metal cations play an important role in the catalytic phenomenon of biomass pyrolysis with three main steps: (1) biomass devolatilization; (2) tar and char formation and (3) tar dissociation by molten salt [START_REF] Ratchahat | Combined molten salt-Ni/Al2O3 as synergistic medium for high-quality syngas production[END_REF].

Several researchers have performed laboratory bench scale experiments using molten salt for biomass pyrolysis. The temperature, amount and composition of molten salt have been shown to influence the pyrolysis reactions. The heating rates of sapwood pyrolysis at 600e900 C in molten salt were 4e10 times higher than those without molten salt, the heating rate can be up to 4100e16700 C/min in a mixture of sodium chloride (NaCl) and potassium chloride (KCl) at 900 C [START_REF] Yasunishi | Wood pyrolysis in molten-salt[END_REF]. Maximum heating rates for cylindrical beech wood pellet (30 mm length and 3.5 mm diameter) pyrolysis in LiNaKF salt were 46e56 C/s [START_REF] Nygard | Thermal history of wood particles in molten salt pyrolysis[END_REF]. The total energy input required for wood pyrolysis at 900 C equal about 20% of its heating value [START_REF] Caubet | Thermochemical conversion of biomass: gasification by flash pyrolysis study[END_REF].

Besides, molten salt changed the proportion of three-state products. When wood pyrolysis in molten carbonates, the yield of liquid products decreased while the yields of gas and char increased [START_REF] Tada | Wood pyrolysis with molten-salt as heating medium[END_REF]. Higher gas yield and less complex bio-oil composition were obtained in a bath of nitrate molten salts (60 wt% NaNO 3 and 40 wt% KNO 3 ) at 500 C [START_REF] Serrano | Thermochemical conversion of C. cardunculus L. in nitrate molten salts[END_REF]. The yields of bio-oil obtained from cellulose/rice stalk pyrolysis in mixtures of metal chlorides and nitrates molten salts at temperatures ranges of 400e600 C increased and then decreased with increasing temperature, and the optimal temperature was 530 C [START_REF] Jiang | Experimental study on thermal pyrolysis of biomass in molten salt media[END_REF]. Lignin was pyrolyzed in mixtures of molten salt (ZnCl 2 and KCl) with molar ratios of 3/7 and 7/6 at temperatures from 500 to 800 C [START_REF] Sada | Pyrolysis of lignins in molten salt media[END_REF]. It was found that a molar ratio of 7/6 at 600 C had a high selectivity to cresols with a maximum yield of 4.6 wt%. The effect of salt amount was also investigated with different salt-to-lignin ratios (1, 2 and 3) [START_REF] Kudsy | Pyrolysis of kraft lignin in the presence of molten ZnCl2-KCl mixture[END_REF]. In addition, the lowest yields of gas products and the highest yield of phenolic compounds except for phenol were obtained at salt-to-lignin ratio of 1. The presence of molten salt in cellulose pyrolysis could increase its conversion to syngas by 30% at 850 C and the maximum conversion of 98 wt% was achieved at 915 C [START_REF] Adinberg | Solar gasification of biomass: a molten salt pyrolysis study[END_REF]. The use of molten salt increased the total hydrocarbon compounds with decreasing oxygenated compounds during coal pyrolysis [START_REF] Rizkiana | Oil production from mild pyrolysis of low-rank coal in molten salts media[END_REF].

Most of these studies were focused on the distribution and composition of products obtained from pyrolysis in molten salt. However, little information is available on the formation mechanism and mutual transformation process of products during solar pyrolysis of biomass in molten salt, despite this information is very important for realizing the high-valued conversion of biomass. Therefore, it is necessary to deeply investigate all products (bio-oil, bio-char and bio-gas) for revealing pathways of intermediates degradation with molten salt. In this study, solar pyrolysis of cotton stalk in molten salt (Li 2 CO 3 -Na 2 CO 3 -K 2 CO 3 ) was conducted in a 4 kW solar simulator reactor. The effects of temperature (450, 550, 650, 750 and 850 C) and molten salt amount (the mass ratio of molten salt to biomass was 0.5, 1, 5 and 10) on the molten salt pyrolysis of biomass were investigated by characterizing the total products.

Experimental

Experimental sample 2.1.1. Cotton stalk

Cotton stalk (CS) was crushed into small particles with size of 0.1e0.4 mm, and then dried at 105 C for 24 h. The proximate analysis of the air-dried cotton stalk feedstock was 77.49 wt% volatile, 4.60 wt% ash, 16.04 wt% fixed carbon, and 1.86 wt% moisture. The elemental composition of the dried feedstock was 46.56 wt% carbon, 6.04 wt% hydrogen, 41.82 wt% oxygen, and 0.79 wt% nitrogen.

Molten salt

Lithium, sodium and potassium carbonates (analytical purity > 99%, Sinopharm Chemical Reagent Co., Ltd, Shanghai, China) were well-mixed with mass ratio of 1:1:1. Then, the ternary eutectic blend salts were put in a corundum crucible and heated to 600 C for 2 h in muffle furnace. The salts were crushed into fine particles by a grinder once they were cooled down to atmospheric temperature. The fine particle salts were defined as molten salts (MS), and the melting point of MS was approximately 400 C [START_REF] Olivares | The thermal stability of molten lithiumesodiumepotassium carbonate and the influence of additives on the melting point[END_REF].

Pyrolysis feedstock

For each pyrolysis experiment, about 1 g of cotton stalk was weighted then well mixed with MS. The mass ratio of MS to cotton stalk for mixture were 0, 0.5, 1, 5, 10 denoted as CS, CS1MS0.5, CS1MS1, CS1MS5 and CS1MS10, respectively.

Experimental procedure 2.2.1. Experimental setup

Solar pyrolysis experiments were performed in molten salt reactor heated by a 4 kW solar simulator instead of the real sun as shown in Fig. 1. Solar simulators can overcome the shortcomings of real sunlight such as the periodicity and randomness of solar irradiance [START_REF] Xiao | Design and characterization of a highflux non-coaxial concentrating solar simulator[END_REF]. What's more, they can present significant advantages in terms of size, cost and operational flexibility [START_REF] Gallo | High flux solar simulators for concentrated solar thermal research: a review[END_REF]. They are very important and valuable for fast and cost-effective indoor laboratory in certain research fields, such as the pyrolysis [START_REF] Arribas | Solar-driven pyrolysis and gasification of low-grade carbonaceous materials[END_REF] and gasification [START_REF] Müller | A pressurized high-flux solar reactor for the efficient thermochemical gasification of carbonaceous feedstock[END_REF] of carbonaceous feedstock.

The solar simulator consists of three 5 kW xenon arc lamps and three accompanying parabolic concentrators. The lamps and concentrator are adjusted to the right position making sure the focus in the same point. The position of solar simulator focus can be changed by adjustable screw devices. The solar simulator can produce concentrated radiation similar to the solar spectral pattern (mostly in the visible and IR spectra) [START_REF] Petrasch | A novel 50 kW 11,000 suns high-flux solar simulator based on an array of xenon arc lamps[END_REF], with irradiative flux densities ranging from 0.2 to 1.5 MW/m 2 within 800 mm diameter circle. A shutter with moving parallel carbon composite blades modulates the concentrated flux impinging molten salt reactor and finally the pyrolysis sample temperature. The molten salt reactor is a 30 mm diameter by 40 mm long cylindrical graphite crucible within a 150 mm transparent Pyrex balloon enclosure, which is positioned vertically at the solar simulator focus. Temperature sensor (K-type thermocouple) is sited at the rear of the reactor in contact with the reactor inner wall to measure the sample temperature. The required heating rate and temperature are set on a PID controller, which controls the shutter opening based on the measured sample temperature. The balloon is swept with nitrogen flow controlled by a mass flowmeter for keeping wall clean.

Experimental procedure

The pyrolysis sample was weighed and placed into the graphite crucible (molten salt reactor). The graphite crucible was put inside the transparent Pyrex balloon and set at the solar simulator focus. The balloon was purged in 3 min by nitrogen gas with a flow rate of 2000 ml/min for inert atmosphere. Then the pyrolysis unit swept by 500 ml/min nitrogen was heated to target temperatures (450, 550, 650, 750 and 850 C) at heating rate of 10 C/s with PID controller on. The pyrolysis temperature was remained stable for 30 min until the reaction was completed. The gas mixture produced from pyrolysis was cooled by the ice-water mixture condensing unit for the condensable volatiles collection, whereas the noncondensing gas was collected with gas bag. After each run, the solid residues were retained in graphite crucible after turning off the solar simulator. The carrier gas flow was still on for 5 min until solid residues cooled to the ambient temperature. Each experiment was run at least three times to ensure the repeatability. The RSD for the experiment was less than 5%. The error bars indicate 95% confidence intervals.

Product recovery and characterization

According to molten salt mass ratio in pyrolysis feedstock, the solid residues was washed with sufficient amount of hydrochloric acid solution (1 mol/L) for removing molten salt. The mixture was stirred in the magnetic agitator for 3 h at room temperature. After stirring, suction filtration device was utilized to filter and wash suspension mixture until the filtrate was neutral. The solid left after filtration was dried at 105 C in oven to constant weight, and then was weighed and recorded as the "bio-char" yield. The "gas" yield was calculated by combining the volumes of all the gases collected in the gas bag. The "bio-oil" yield was calculated by the mass difference of the conical condenser after the reaction.

The bio-char was weighed and characterized as our previous study [START_REF] Zeng | Effects of temperature and Mg-based additives on properties of cotton stalk torrefaction products[END_REF]. Ultimate analysis was carried out in a CHNS/O elementary analyzer (Vario Micro Cube, Germany). The specific surface area and pore size distribution were analyzed with Surface Area and Porosity analyzer (TriStar II, Micromeritics Instrument Corp., Norcross, USA). Prior to each analysis, the sample was degassed at the temperature of 150 C under a vacuum for 20 h. Nitrogen adsorption isotherm was obtained at the temperature of À196 C in the relative pressures (P/P 0 ) range between 0.05 and 0.99. The BET surface area, micropore volume and total pore volume were calculated with the BrunauereEmmetteTeller (BET) equation, t-plot method and single point adsorption total pore volume analysis, respectively [START_REF] Chen | Biomass-based pyrolytic polygeneration system on cotton stalk pyrolysis: influence of temperature[END_REF].

The gas compositions were analyzed using a dual-channel gas chromatography (Panna A91, Panna Instruments Co., Ltd, Changzhou, China). H 2 , N 2 , O 2 , CO 2 , CH 4 and CO were separated by Porapak Q column and 5 A zeolite molecular sieve column with He as carrier gas and analyzed by a thermal conductivity detector (TCD). Olefins and alkanes (C 2þ ) were separated using a Al 2 O 3 /KCl capillary column with N 2 as carrier gas and analyzed by a flame ionization detector (FID). Content of N 2 in produced gas was calculated as the reference to obtain the yield of each gas. Heating program of columns was set as below, columns was maintained at 80 C for 4 min and heated to 150 C at a heating rate of 10 C/min, then maintained at 150 C for 6 min.

The main compositions of bio-oil were determined using gas chromatography-mass spectrometry (GC-MS; Agilent 7890 A series GC coupled with an HP5975 MS detector) with He. The detail of GC-MS and test method has already been described in our previous study [START_REF] Zeng | Influence of torrefaction with Mg-based additives on the pyrolysis of cotton stalk[END_REF]. All the compounds were identified based on NIST mass spectral data library with response factors for GC quantification determined by internal standard method [START_REF] Zeng | Characterization of solar fuels obtained from beech wood solar pyrolysis[END_REF].

Results and discussion

Final product distribution

Fig. 2 compares the products distribution from solar pyrolysis of cotton stalk without and with molten salt. For CS pyrolysis (without molten salt), the gas yield noticeably increased from 24.09 wt% to 41.35 wt% when pyrolysis temperature increased from 450 to 850 C shown in Fig. 2a, while the liquid yield significantly decreased from 42.08 wt% to 33.65 wt% indicated in Fig. 2b, char yield decreased from 34.64 wt% to 25.11 wt% depicted in Fig. 2c. This result suggests that the increase in the gas yield is primarily ascribed to the secondary tar cracking and char decomposition, which can be confirmed by the reduce in liquid and char yields [START_REF] Zeng | Solar pyrolysis of beech wood: effects of pyrolysis parameters on the product distribution and gas product composition[END_REF]. The secondary reactions of tar into gas prevail at temperatures above 600 C [START_REF] Neves | Characterization and prediction of biomass pyrolysis products[END_REF]. The changing trend of product distribution with pyrolysis temperature was almost the same for cotton stalk pyrolysis with and without molten salt. For CS1MS10 pyrolysis, the gas yield distinctly increased from 31.07 wt% to 82.57 wt% with pyrolysis temperature rising from 450 to 850 C (Fig. 2a), while the liquid yield decreased from 35.67 wt% to 28.68 wt% (Fig. 2b), char yield remarkably decreased from 38.50 wt% to 19.50 wt% (Fig. 2c).

For the same pyrolysis temperature, the use of molten salt as pyrolysis media increased gas yield at the expense of liquid product or char. The increase of gas yield was strongly due to the catalytic activity of the alkali species in the molten salt for cotton stalk pyrolysis. Molten salt made sure a fast heat transfer rate because of its high thermal conductivity and heat capacity [START_REF] Hathaway | Solar gasification of biomass: kinetics of pyrolysis and steam gasification in molten salt[END_REF]. Then the cotton stalk thorough dispersion in liquid salt enabled rapid heating for higher pyrolysis rate [START_REF] Adinberg | Solar gasification of biomass: a molten salt pyrolysis study[END_REF]. It was expected to obtain more primary vapor products due to more rapid devolatilization at pyrolysis [START_REF] Zeng | The effect of temperature and heating rate on char properties obtained from solar pyrolysis of beech wood[END_REF]. With the help of alkali carbonate salts catalysis effect, the primary vapor secondary cracking reaction into gas product would be promoted at higher pyrolysis temperature [START_REF] Nzihou | A review of catalysts for the gasification of biomass char, with some reference to coal[END_REF]. Meanwhile, the char and molten slat might react to produce more gas at higher temperature (>700 C) [START_REF] Lee | Gasification of bamboo carbon with molten alkali carbonates[END_REF], and this also made the mass balance of pyrolysis over 100% at high temperature (750 and 850 C). Then the increase of gas yield became more significant with increasing pyrolysis temperature indicated in Fig. 2a.

At pyrolysis temperature of 850 C, the gas yield continued to rise from 64.31 wt% to 82.57 wt% when mass ratio of molten salt to cotton stalk increased from 0.5 to 10 shown in Fig. 2a. A higher molten salt mass ratio meant more effective contact of alkali carbonates with the cotton stalk so that the secondary reaction of tar into gas product increased [START_REF] Rizkiana | Oil production from mild pyrolysis of low-rank coal in molten salts media[END_REF]. Meanwhile, more molten salt could participate in the reaction with char, thus producing more gas. However, there was no statistically significant change in gas yield with different molten salt mass ratios when the pyrolysis temperatures were lower than 850 C. It was mainly attributed to the not so obvious catalytic effect of molten salt (Na 2 CO 3 and K 2 CO 3 ) for tar secondary reactions and the reaction with char when the temperature was lower than 800 C [START_REF] Adinberg | Solar gasification of biomass: a molten salt pyrolysis study[END_REF].

The existence of molten salt could not only catalyze the decomposition of volatile matter, but also might inhibit the release of volatile matter, thus reducing the yield of the bio-oil and increased the yield of bio-char. Yasunishi, A. et al. [START_REF] Yasunishi | Wood pyrolysis in molten-salt[END_REF] has also found that when wood was pyrolyzed in molten carbonates, the yield of liquid products decreased while the yields of gas and char increased. The inhibitory effect of molten salt on volatile release enhanced with the increase of molten salt mass ratio. So the liquid yields obtained from CS1MS5 and CS1MS10 were lower than those from CS1MS0.5 and CS1MS1 indicated in Fig. 2b.

It was observed that molten salt use slightly increased the char yield at low pyrolysis temperatures (450e650 C) but it further decreased when temperature increased to 850 C (Fig. 2c). The displacement of pyrolysis reaction at lower temperatures induced by alkali metals favors the charring reactors versus fragmentation and depolymerization reactions in the primary decomposition of biomass, thereby increasing char yield [START_REF] Blasi | Influences of the chemical state of alkaline compounds and the nature of AlkaliMetal on wood pyrolysis[END_REF]. The decrease of char yield at higher temperatures (750 and 850 C) were mainly because K 2 CO 3 and Na 2 CO 3 could further promote the cleavage of C-H and C-O bonds and aromatization of char during pyrolysis [START_REF] Hu | Effects of inherent alkali and alkaline earth metallic species on biomass pyrolysis at different temperatures[END_REF] while carbon in the char would simultaneously reacted with carbonates as the reaction (1) [START_REF] Nagase | Kinetics and mechanisms of the reverse Boudouard reaction over metal carbonates in connection with the reactions of solid carbon with the metal carbonates[END_REF]. However, considering the error bars, there was no big difference of the char yields obtained with different molten salt mass ratios under the same pyrolysis temperature.

CO 2À 3 þ C/2CO þ O 2À (2

Pyrolysis gas composition

Fig. 3 shows the yields of H 2 , CO 2 , CO, CH 4 and C 2þ during pyrolysis of cotton stalk without and with molten salt at different temperatures. For CS pyrolysis (without molten salt), the H 2 releasing clearly increased from 2.15 to 81.64 ml/g with pyrolysis temperature increasing from 450 to 850 C shown in Fig. 3a. And CO distinctly increased from 42.07 to 138.55 ml/g depicted in Fig. 3b, CH 4 noticeably increased from 8.93 to 47.86 ml/g indicated in Fig. 3c, C 2þ increased from 5.00 to 18.68 ml/g as shown in Fig. 3d. While, temperature had a limited influence on the yield of CO 2 in this temperature range (Fig. 3e). Below 450e550 C, most of permanent gases directly came from degradation of biomass chemical structures (hemicellulose, cellulose and lignin) as volatiles secondary reactions was negligible [START_REF] Neves | Characterization and prediction of biomass pyrolysis products[END_REF]. However, the increase of H 2 and CO were mainly due to tar secondary reactions as temperature increases above 450e550 C [START_REF] Morf | Mechanisms and kinetics of homogeneous secondary reactions of tar from continuous pyrolysis of wood chips[END_REF]. Simultaneously, most of the remaining tars convert into CH 4 and C 2þ [START_REF] Zeng | Characterization of solar fuels obtained from beech wood solar pyrolysis[END_REF].

The changing trend of gas compositions with pyrolysis temperature was almost the same for cotton stalk pyrolysis with and without molten salt. For CS1MS10 pyrolysis, H 2 yield clearly increased from 11.33 to 162.43 ml/g with pyrolysis temperature increasing from 450 to 850 C shown in Fig. 3a. And CO distinctly increased from 36.60 to 388.68 ml/g depicted in Fig. 3b noticeably increased from 13.96 to 43.07 ml/g indicated in Fig. 3c, C 2þ increased from 8.30 to 18.81 ml/g as observed in Fig. 3d. At the same pyrolysis temperature, H 2 yield showed a general increasing trend with molten salt media (Fig. 3a). With Na containing molten salt, the added Na þ and eCOOH or eOH groups in cotton stalk reacted as reactions ( 2) and ( 3) for releasing H þ , which promoted H 2 formation [START_REF] Ding | Catalytic effects of Na2CO3 additive on coal pyrolysis and gasification[END_REF]. Furthermore, the increase trend of H 2 yield became more pronounced with pyrolysis temperature rising. It was mainly due to the enhancing catalytic activity of the alkali species in the molten salt for tar cracking reaction (4) at higher temperatures [START_REF] Hathaway | Solar gasification of biomass: kinetics of pyrolysis and steam gasification in molten salt[END_REF].

Na þ þ À COOH/ À COONa þ H þ (3) Na þ þ À OH/ À ΟNa þ H þ (4)  C x H y O z à / C þ Tar þ CO þ H 2 þ CH 4 þ C n H m ( 5
)
Using molten salt as pyrolysis media also increased CO yield by 1e3 times compared to those without molten salt when temperature was above 550 C. The catalytic cracking of tar at higher temperatures for producing CO as (4) was enhanced by molten salt [START_REF] Ding | Catalytic effects of Na2CO3 additive on coal pyrolysis and gasification[END_REF]. With the use of molten salt (Na 2 CO 3 and K 2 CO 3 ), CO 2 produced during primary pyrolysis stage would react with carbon bonded with Na or K as reactions ( 5), ( 6) and ( 7) before its evolution from biomass or coke surface [START_REF] Rizkiana | Oil production from mild pyrolysis of low-rank coal in molten salts media[END_REF]. What's more, the char would react with carbonates over 700 C for generating CO by an intermolecular redox reaction (1) [START_REF] Nagase | Kinetics and mechanisms of the reverse Boudouard reaction over metal carbonates in connection with the reactions of solid carbon with the metal carbonates[END_REF]. The extend of these reactions increased with temperature rising, and therefore CO yield increase became more and more significant with increasing temperature.

M 2 CO 3 þ 2C/2M þ 3CO (6) 2M þ CO 2 /M 2 O þ CO (7) M 2 O þ CO 2 /M 2 CO 3 ( 8
)
In contrast, a lower CH 4 production was in the molten salt pyrolysis when temperature was above 550 C (Fig. 3c). It was consistent with a report that 30% of CH 4 conversion at 900 C in molten salt as [START_REF] Zeng | Modeling of beech wood pellet pyrolysis under concentrated solar radiation[END_REF] [START_REF] Xianquan | Reduction behavior of methane in alkali molten carbonates[END_REF]. A comparison between pyrolysis without and with molten salt showed that variation on the yields of C 2þ were very small (Fig. 3d). It was also reported by Blasi et al. [START_REF] Blasi | Influences of the chemical state of alkaline compounds and the nature of AlkaliMetal on wood pyrolysis[END_REF]. Significant increase was observed for CO 2 with the use of molten salt (Fig. 3e). The increase of CO 2 was mainly due to CO conversion via water gas shift reaction [START_REF] Zeng | Combined effects of initial water content and heating parameters on solar pyrolysis of beech wood[END_REF] and partly because of carbonate salt decomposition reactions [START_REF] Nygård | Review of thermal processing of biomass and waste in molten salts for production of renewable fuels and chemicals[END_REF], [START_REF] Hathaway | Solar gasification of biomass: kinetics of pyrolysis and steam gasification in molten salt[END_REF] [START_REF] Hathaway | Steam gasification of plant biomass using molten carbonate salts[END_REF][START_REF] Rizkiana | Oil production from mild pyrolysis of low-rank coal in molten salts media[END_REF].

CH 4 / C þ 2H 2 (9) CO þ H 2 O/CO 2 þ H 2 (10) M 2 CO 3 þ H 2 O/2MOH þ CO 2 (11) M 2 CO 3 / M 2 O þ CO 2 ( 12
)
It indicated that the increase of mass ratio of molten salt to cotton stalk from 0.5 to 10 could promote H 2 and CO production at pyrolysis temperature of 850 C. Increasing molten salt content increased the secondary reaction of tar since more tar could contact with the salt. However, there was almost no change in H 2 and CO yields with different molten salt mass ratios when the pyrolysis temperatures were lower than 850 C. As shown in Fig. 3c ande, the yields of CH 4 and CO 2 obtained with 10 and 5 mass ratio of molten salt to cotton stalk was higher than those with 0.5 and 1 mass ratio of molten salt to cotton stalk. This may be because more alkali metal salts (Na, K) could promote the cracking and catalytic reforming reactions between organic compounds and methane [START_REF] Zhou | Hydropyrolsis characteristics and kinetics of potassium-impregnated pine wood[END_REF].

Bio-oil characterization

Fig. 4 describes the composition of bio-oil from cotton stalk pyrolysis without and with molten salt at temperatures of 450e850 C. Bio-oil contained a wide range of organic compounds that can be divided into several main groups: Acids, Linear alcohols/ aldehydes/ketones, Furans, Cyclopentanones, Phenols, Aromatics and N-species. For cotton stalk pyrolysis (Fig. 4a), the relative content of aromatics tended to increase from 10.07% to 45.30% with the pyrolysis temperature from 450 to 850 C. Significantly higher aromatic hydrocarbons could be attributed to the enhanced formation of polycyclic aromatic hydrocarbons (PAHs) such as phenanthrene, anthracene, naphthalene with temperature rising [START_REF] Akhtar | A review on operating parameters for optimum liquid oil yield in biomass pyrolysis[END_REF]. While acids, linear alcohols/aldehydes/ketones and phenols gradually reduced from temperature 450e850 C. This was because these functional groups were mainly formed from the depolymerization of hemicellulose and cellulose at temperature lower than 350 C [START_REF] Collard | A review on pyrolysis of biomass constituents: mechanisms and composition of the products obtained from the conversion of cellulose, hemicelluloses and lignin[END_REF]. As the final pyrolysis temperature increased, the functional groups experienced secondary cracking for producing more stable species.

For cotton stalk pyrolysis with molten salt (Fig. 4b), the changing trend of bio-oil composition with pyrolysis temperature was very different from that without molten salt. Molten salt pyrolysis caused a noticeable decrease in bio-oil acids content (3.14%e0.77% versus 31.67%e22.26% of original cotton stalk) among pyrolysis temperature range of 450e850 C. The decrease of acids was mainly due to enhanced catalytic cracking of the cracking of carboxylic acids by molten salt components (Li 2 CO 3 -Na 2 CO 3 -K 2 CO 3 ) [START_REF] Wang | Catalytic conversion of carboxylic acids in bio-oil for liquid hydrocarbons production[END_REF]. Moreover, alkali and alkaline earth metals in molten salt were able to promote the catalytic conversion of acidic substances and oxygenates in the bio-oil through ketonization and the aldol reaction for more ketones production [START_REF] Stefanidis | In-situ upgrading of biomass pyrolysis vapors: catalyst screening on a fixed bed reactor[END_REF]. However, the cyclopentanones content in molten salt pyrolysis was more than twice that of pyrolysis without molten salt. The addition of Li, Na or K in molten salt would be attributed to the improvement of cyclopentanones formation, as they would enhance the cleavage of C-O bond and the formation of new C-C bond in bio-oil [START_REF] Zhou | Selective conversion of furfural to cyclopentanone over CNT-supported Cu based catalysts: model reaction for upgrading of bio-oil[END_REF]. When cotton stalk was pyrolyzed in molten salt, phenols significantly decreased. In the presence of alkali and alkaline earth metal catalyst, lignin-derived phenolics undergo oxygen-aromatic carbon bond cleavage to form aromatic hydrocarbons [START_REF] Mourant | Mallee wood fast pyrolysis: effects of alkali and alkaline earth metallic species on the yield and composition of bio-oil[END_REF][START_REF] Zabeti | In situ catalytic pyrolysis of lignocellulose using alkali-modified amorphous silica alumina[END_REF]. It was also confirmed by the fact that molten salt presence increased the relative content of aromatics in the bio-oil.

The main functional groups (Acids, O-species (furans/aldehydes/ ketones/alcohols), Phenols and Aromatics) are selected for further comparison indicated in Fig. 5. As depicted in Fig. 5a, the reduction of acids in pyrolysis bio-oil was most evident when the mass ratio of molten salt to cotton stalk was 10. The most obvious decrease of acid was mainly due to optimum catalytic effect achieved with this mass ratio [START_REF] Jiang | Experimental study on thermal pyrolysis of biomass in molten salt media[END_REF]. However, with the mass ratio of molten salt to cotton stalk increased from 0.5 to 10, the relative content of acids in pyrolysis bio-oil gradually decreased. As shown in Fig. 5b, the ospecies (furans/aldehydes/ketones/alcohols) increased when mass ratio of molten salt to cotton stalk increased from 0.5 to 10 at lower pyrolysis temperature (450, 550 and 650 C). Ketonization and the aldol reactions for more ketones and alcohols production were promoted due to the effective contact of bio-oil compounds with more molten salt catalyst at high mass ratio of molten salt to cotton stalk. However, the o-species might be transformed into gaseous products, which would reduce its relative content at high temperature. Increasing molten salt mass ratio promoted the catalytic cracking of phenols when pyrolysis temperature above 550 C indicated in Fig. 5c. It should be noted from Fig. 5d that there was positive correlation between the increased content of aromatics and mass ratio of molten salt to cotton stalk (from 0.5 to 10). In fact, the higher the molten salt ratio was used, the better dispersion of cotton stalk in molten salt obtained. Then increasing molten salt content increased the secondary tar reaction of into aromatics since it enhanced the catalytic effect of molten salt (Na 2 CO 3 and K 2 CO 3 ) [START_REF] Rizkiana | Oil production from mild pyrolysis of low-rank coal in molten salts media[END_REF].

Bio-char characterization

The elemental analysis of bio-char obtained from pyrolysis of cotton stalk without and with molten salt are listed in Table 1. For cotton stalk pyrolysis, the carbon content in bio-char increased with increasing pyrolysis temperature while the hydrogen and oxygen contents declined sharply. The carbonaceous structure graphitization at high temperatures lead to increased the carbon content [START_REF] Zeng | The effect of temperature and heating rate on char properties obtained from solar pyrolysis of beech wood[END_REF]. The decrease of hydrogen and oxygen contents was due to the char weak bonds cracking and cleavage [START_REF] Kim | Influence of pyrolysis temperature on physicochemical properties of biochar obtained from the fast pyrolysis of pitch pine (Pinus rigida)[END_REF]. The changing trend of bio-char composition with pyrolysis temperature was different for cotton stalk pyrolysis with and without molten salt. At the same pyrolysis temperature, carbon content showed a general decreasing trend while oxygen increased with molten salt media. The decrease of carbon content was mainly due to the enhanced C-C bonds cleavage by alkaline metals in molten salt [START_REF] Blasi | Influences of the chemical state of alkaline compounds and the nature of AlkaliMetal on wood pyrolysis[END_REF]. Alkaline cations due to small size could penetrate into the biomass/ char structure and cause the bond of oxygen as reactions ( 2) and ( 3) [START_REF] Ding | Catalytic effects of Na2CO3 additive on coal pyrolysis and gasification[END_REF]. However, there was no obvious influence of molten salt mass ratio on bio-char composition. Table 2 lists the BET surface area and pore volume of bio-char obtained from pyrolysis of cotton stalk without and with molten salt. For cotton stalk pyrolysis, the BET surface area increased from 2.12 to 222.31 m 2 /g as temperature increased from 450 to 850 C. Meanwhile, the micropore volume increased from 0.0011 to 0.0978 cm 3 /g, the total pore volume increased from 0.0038 to 0.1395 cm 3 /g. Intensifying volatile release with increasing temperature caused the increase of surface area and pore volume [START_REF] Kim | Influence of pyrolysis temperature on physicochemical properties of biochar obtained from the fast pyrolysis of pitch pine (Pinus rigida)[END_REF]. The changing trend of bio-char BET surface area and pore volume with pyrolysis temperature was almost the same for cotton stalk pyrolysis with and without molten salt. At the same pyrolysis temperature, BET surface and pore volume showed an obvious increasing trend with molten salt media. For CS1MS10 pyrolysis, the BET surface area increased from 1.90 m 2 /g to 873.44 m 2 /g as temperature increased from 450 to 850 C. And the micropore volume increased from 0.0009 to 0.3210 cm 3 /g, the total pore volume increased from 0.0051 to 0.5567 cm 3 /g. At 850 C, the BET surface area of CS1MS10 bio-char increased about 4 times compared with that of CS bio-char. Simultaneously, both of micropore and total pore volumes increased about 4 times. The use of molten salt catalyzed the pyrolysis reaction of cotton stalk, which enhanced volatile release for creating bio-char with higher BET surface area and pore volume. In addition, the generated H 2 O and CO 2 during pyrolysis process reacted with carbon for porous carbon formation, which called activation process [START_REF] Yin | Harvesting capacitive carbon by carbonization of waste biomass in molten salts[END_REF]. Molten salt just could promote the activation process for more porous bio-char with high surface area and pore volume. However, bio-char obtained from C1MS5 pyrolysis at 850 C had the highest surface area of 972.57 m 2 /g and the biggest total pore volume of 0.6203 cm 3 /g. BET surface area and pore volume of bio-char decreased with mass ratio of molten salt to cotton stalk increasing from 5 to 10. The results proved that loading excessive amounts of molten salt could enter the inner interspace through the surface crack of char, which resulted in blocking some micropores and decreasing the specific surface of bio-char [START_REF] Ding | Catalytic effects of Na2CO3 additive on coal pyrolysis and gasification[END_REF].

Conclusion

Solar pyrolysis of cotton stalk in molten salt media has been investigated with a 4 kW solar simulator. At pyrolysis temperature of 850 C, the syngas mainly composed of CO and H 2 , continued to rise from 41.35 wt% to 82.57 wt% when mass ratio of molten salt to cotton stalk increased from 0 to 10. Pyrolysis in molten salt significantly decreased bio-oil acids (3.14%e0.77% versus 31.67%e 22.26% of original cotton stalk) and phenols, while increased aromatics among pyrolysis temperature range of 450e850 C. There was positive correlation between the increased content of aromatics and mass ratio of molten salt to cotton stalk (from 0.5 to 10). The bio-char carbon content showed a general decreasing trend while oxygen, BET surface and pore volume increased with using molten salt as pyrolysis media. BET surface area and pore volume of bio-char decreased with mass ratio of molten salt to cotton stalk increasing from 5 to 10. Solar pyrolysis of biomass in molten salt can be considered as a promising method for producing syngas (rich in H 2 and CO), high quality bio-oil (little acids) and bio-char (large BET surface area), which storing solar energy in chemical form at the same time.

Fig. 1 .

 1 Fig. 1. Schematic of the solar pyrolysis in molten salt experimental setup.

Fig. 2 .

 2 Fig. 2. Effects of temperature and molten salt amount on the product distribution of solar pyrolysis: (a) Bio-gas yield, (b) Bio-oil yield, and (c) Bio-char yield.

Fig. 3 .

 3 Fig. 3. Effects of temperature and molten salt amount on the gas product compositions of solar pyrolysis: (a) H 2 yield, (b) CO yield, (c) CH 4 yield, (d) C 2þ yield, and (e) CO 2 yield.

Fig. 4 .

 4 Fig. 4. Effects of temperature on the bio-oil chemical composition from solar pyrolysis of: (a) CS and (b) CS1MS1.

Fig. 5 .

 5 Fig. 5. Effects of temperature and molten salt amount on the main functional group distribution of bio-oil from solar pyrolysis: (a) Acetic acid, (b) O-species (furans/aldehydes/ ketones/alcohols), (c) Phenols, and (d) Aromatics.

Table 1

 1 Effect of pyrolysis temperature and mass ratio of molten salt to biomass on char composition.

	Char samples	Ultimate analysis (wt.%, daf)
		C	H	N	O a
	CS-450	81.97	3.51	1.89	12.63
	CS-550	86.06	3.10	1.76	9.08
	CS-650	90.13	2.45	1.44	5.98
	CS-750	93.42	1.96	1.31	3.31
	CS-850	94.42	1.56	1.34	2.69
	CS1MS0.5-550	83.70	2.66	1.27	12.37
	CS1MS0.5-850	84.92	1.02	0.62	13.44
	CS1MS1-550	75.79	2.59	1.19	20.43
	CS1MS1-850	81.85	1.75	0.28	16.12
	CS1MS5-550	77.86	3.60	1.27	17.27
	CS1MS5-850	83.71	1.25	0.16	14.88
	CS1MS10-450	74.16	3.51	1.29	21.04
	CS1MS10-550	78.54	2.61	1.23	17.63
	CS1MS10-650	78.44	2.31	1.01	18.24
	CS1MS10-750	79.16	2.11	0.45	18.28
	CS1MS10-850	78.98	1.53	0.09	19.40
	daf: dry ash free basis.

a By difference.

Table 2

 2 Effect of pyrolysis temperature and mass ratio of molten salt to biomass on char surface area and pore volume.

	Char samples	BET surface area (m 2 /g)	Pore volume of micropore (cm 3 /g)	Total pore volume (cm 3 /g)
	CS-450	2.12	0.0011	0.0038
	CS-550	11.63	0.0155	0.0155
	CS-650	84.24	0.0320	0.0556
	CS-750	192.26	0.0854	0.1157
	CS-850	222.31	0.0978	0.1395
	CS1MS0.5-550	5.68	0.0040	0.0087
	CS1MS0.5-850	802.29	0.2935	0.4900
	CS1MS1-550	1.41	0.0018	0.0041
	CS1MS1-850	904.60	0.3128	0.5387
	CS1MS5-550	1.94	0.0027	0.0070
	CS1MS5-850	972.57	0.2881	0.6203
	CS1MS10-450	1.90	0.0009	0.0051
	CS1MS10-550	6.08	0.0035	0.0080
	CS1MS10-650	340.41	0.1533	0.2060
	CS1MS10-750	641.49	0.2784	0.3895
	CS1MS10-850	873.44	0.3210	0.5567

Acknowledgement

This work was supported by the National Natural Science Foundation of China (Project No. 51706083 and 51622604), and Natural Science Foundation of Shenzhen (JCYJ20170818164006890). We also express appreciation of the technical support from Analytical and Testing Center in Huazhong University of Science & Technology (http://atc.hust.edu.cn).


Gaétan Bardy
email: bardy@gmail.com

Pierre Biver
email: pierre.biver@total.com

Guillaume Caumon
email: guillaume.caumon@univ-lorraine.fr

Philippe Renard
email: philippe.renard@unine.ch

Oil production uncertainty assessment by predicting reservoir production curves and confidence intervals from arbitrary proxy responses

Keywords: Uncertainty Quantification, Model Ranking, Geostatistics, Ensemble Methods. Highlights

Underground fluid flow in hydrocarbon reservoirs (or aquifers) is difficult to predict accurately due to geological and petrophysical uncertainties. To quantify that uncertainty, several spatial statistical methods are often used to generate an ensemble of subsurface models representing and sampling these uncertainties. However, to predict the uncertainties in terms of flow responses, one needs to run a forward flow simulator (often multiphase flow in transient state) on every model of this ensemble and this generally entails intractable computational costs. Approximate solutions (flow proxies) can help addressing this challenge but introduce physical simplifications whose impact on the uncertainty quantification is difficult to characterize. This paper proposes a workflow to assess the dynamic reservoir behavior uncertainties from an input ensemble of realizations sampling geological and geophysical uncertainties. Analytical reservoir production curves are estimated from proxy distances computed between all ensemble members and from a few accurate flow responses computed on a subset of the ensemble. A randomization process accounting for proxy quality and for model selection is used to assess confidence intervals about reservoir production quantile curves. The process can use both static and dynamic proxies and also permits to compare their efficiency. A case study on a real turbiditic reservoir shows the applicability of the method, and highlights the value of even a simple proxy to increase the confidence about future reservoir production.

Introduction

Dynamic reservoir simulations are used to help decision making in the oil industry. Advances in computational methods make complex flow simulation feasible nowadays on detailed deterministic reservoir models honoring subsurface data (e.g., [START_REF] Obi | Giga Cell Compositional Simulation[END_REF][START_REF] Wang | A scalable parallel black oil simulator on distributed memory parallel computers[END_REF]. However, there is often a practical gap between the number of flow simulations that engineers would need to run to solve a field management problem and the actual number of simulations that can be managed with available computational resources. Indeed, a large number of simulations are needed both to test multiple production scenarios and to account for subsurface uncertainty (see for instance [START_REF] Caers | Modeling Uncertainty in the Earth Sciences[END_REF][START_REF] Corre | Integrated Uncertainty Assessment For Project Evaluation and Risk Analysis[END_REF][START_REF] Fetel | Reservoir flow uncertainty assessment using response surface constrained by secondary information[END_REF][START_REF] Gross | Response surface approaches for large decision trees: decision making under uncertainty[END_REF][START_REF] Manceau | Use of experimental design methodology to make decisions in an uncertain reservoir environment from reservoir uncertainties to economic risk analysis[END_REF][START_REF] Maschio | A new methodology to reduce uncertainties in reservoir simulation models using observed data and sampling techniques[END_REF][START_REF] Schiozer | Use of representative models in the integration of risk analysis and production strategy definition[END_REF][START_REF] Subbey | Prediction under uncertainty in reservoir modeling[END_REF][START_REF] Xie | Integration of Shale-Gas-Production Data and Microseismic for Fracture and Reservoir Properties With the Fast Marching Method[END_REF]. In this paper, we focus on the specific problem of capturing efficiently the uncertainties on the temporal evolution of field production under a fixed reservoir production scenario. The knowledge of these uncertainties may affect development decisions of green fields (reservoirs not yet in production where no historical data are available), and also provides important information for narrowing the search space during production history matching of brown fields (mature fields with several years of production history). As geological and geophysical interpretations may take many forms, we choose to incorporate this expert interpretive knowledge using a large ensemble of detailed static 3-D reservoir models.

Accurate flow simulation on each of these models being computationally prohibitive, one may need to select only a few models from the initial ensemble to make flow forecasts. However, statistics derived only from a few samples are not robust. This calls for carefully selecting a set of representative models from the full ensemble using some prior ranking of models [START_REF] Ballin | Prediction Of Uncertainty In Reservoir Performance Forecast[END_REF][START_REF] Deutsch | Fortran programs for calculating connectivity of three-dimensional numerical models and for ranking multiple realizations[END_REF][START_REF] Yazdi | Fast screening of geostatistical realizations for SAGD reservoir simulation[END_REF][START_REF] Rahim | Reservoir Geological Uncertainty Reduction: an Optimization-Based Method Using Multiple Static Measures[END_REF] and using specific methodologies such as bootstrap to assess confidence intervals around reservoir responses [START_REF] Scheidt | Bootstrap confidence intervals for reservoir model selection techniques[END_REF]. To-date, several approximate solutions (proxies) have been proposed to help quickly select representative models. Upscaling computes equivalent petrophysical properties at a coarser scale than the initial detailed model [START_REF] Durlofsky | Upscaling and gridding of fine scale geological models for flow simulation[END_REF]. Alternatively, simplified physics such as streamline simulation [START_REF] Thiele | Simulating flow in heterogeneous systems using streamtubes and streamlines[END_REF] may be used on the fine-scale reservoir models. The computational cost of upscaling and streamline simulation and the relative complexity of these methods have motivated even simpler approximations of the dynamic reservoir response based on connectivity [START_REF] Alabert | Stochastic Models of Reservoir Heterogeneity: Impact on Connectivity and Average Permeabilities[END_REF][START_REF] Ballin | Prediction Of Uncertainty In Reservoir Performance Forecast[END_REF][START_REF] Deutsch | Fortran programs for calculating connectivity of three-dimensional numerical models and for ranking multiple realizations[END_REF][START_REF] Renard | Connectivity metrics for subsurface flow and transport[END_REF] or Fast Marching methods [START_REF] Hovadik | Predicting waterflood behavior by simulating earth models with no or limited dynamic data: From model ranking to simulating a billion-cell model[END_REF][START_REF] Xie | Integration of Shale-Gas-Production Data and Microseismic for Fracture and Reservoir Properties With the Fast Marching Method[END_REF]. To investigate how far it is possible to go with the simplification (and therefore with the reduction of the computing cost) we can even consider using an extremely simple proxy such as the Stock Tank Original Oil In Place (STOOIP). In this case, we will obtain only a scalar value and will not get an approximate recovery curve from this computation, but we will still be able to use it as a proxy in the proposed methodology.

In all these cases, the surrogate or proxy responses cannot or should not be used directly to characterize production and the corresponding uncertainties. Instead, it should somehow be calibrated against accurate flow simulations on a few models. Indeed, inaccuracies in the proxy results may lead to biased reservoir forecasts and a possible underestimation of the associated uncertainties [START_REF] Doherty | Use of paired simple and complex models to reduce predictive bias and quantify uncertainty[END_REF][START_REF] Josset | Local and Global Error Models to Improve Uncertainty Quantification[END_REF]Josset et al., 2015).

For this reason, combining proxy evaluation on all models and accurate flow simulation on a few models seems to provide a good compromise to obtain representative uncertainty assessment while harnessing computational costs [START_REF] Doherty | Use of paired simple and complex models to reduce predictive bias and quantify uncertainty[END_REF]Effendiev et al., 2009;Ginsgourger et al., 2013;[START_REF] Josset | Local and Global Error Models to Improve Uncertainty Quantification[END_REF]Josset et al., 2015;[START_REF] Maschio | Bayesian history matching using artificial neural network and Markov Chain Monte Carlo[END_REF][START_REF] Scheidt | Representing Spatial Uncertainty Using Distances and Kernels[END_REF][START_REF] Scheidt | A multi-resolution workflow to generate high-resolution models constrained to dynamic data[END_REF]). Yet, the choice of a particular proxy can be consequential in practical case In this work, we build on a class of methods which exploit distances between reservoir models, initially published by Suzuki and Caers (2008). These distances can be defined from many possible proxies, including simple scalar evaluations such as STOOIP, and do not explicitly call for a particular type of parameterization of the static models. For instance, these approaches can be applied to reservoir models of different geometries (Suzuki et al., 2008). The core idea is to define a distance which is correlated to the flow response, then to use this distance to explore the model space efficiently. Based on these concepts, [START_REF] Scheidt | Representing Spatial Uncertainty Using Distances and Kernels[END_REF] proposed a general model selection technique which forms the basis of the present paper and of other recent works dealing with reservoir uncertainty assessment [START_REF] Josset | Local and Global Error Models to Improve Uncertainty Quantification[END_REF]Josset et al., 2015a;[START_REF] Scheidt | Bootstrap confidence intervals for reservoir model selection techniques[END_REF][START_REF] Scheidt | Bayesian Evidential Learning, in: Quantifying Uncertainty in Subsurface Systems[END_REF] and history matching [START_REF] Ginsbourger | Distance-based kriging relying on proxy simulations for inverse conditioning[END_REF]Josset et al., 2015b;[START_REF] Scheidt | A multi-resolution workflow to generate high-resolution models constrained to dynamic data[END_REF]. Essentially, this class of approaches computes distances from proxies between all possible models, then uses Multi-Dimensional Scaling (MDS) to map models in a feature space where clustering and model selection is performed. Accurate flow simulation is then performed on this small subset to assess reservoir performance uncertainty. To account for limited sample size, [START_REF] Scheidt | Bootstrap confidence intervals for reservoir model selection techniques[END_REF] propose a bootstrap method to assess the uncertainties about reservoir responses derived from these small set of samples. [START_REF] Josset | Local and Global Error Models to Improve Uncertainty Quantification[END_REF] extend this methodology to estimate true responses for all the models from multi-scale flow simulation on all the models in order to obtain more robust uncertainty assessments and confidence intervals about tracer breakthrough curves in aquifers. Their method is partly similar to ours, but only considers flow simulation proxies which account for time and space explicitly, as it relies on the comparison of approximate and exact breakthrough curves at each simulation time step. Our method, in contrast, can use the simpler and faster proxies discussed above such as connectivity-based proxies and static proxies such as STOOIP.

Our main contribution is to propose a methodology applicable to all classes of reservoir proxies. For this, we present a technique able to generate surrogate analytical production curves for all models of the initial ensemble based on the full set of proxy distances and a few accurate flow simulations on a few reference models. After defining the problem in more formal terms, this paper presents an overview of the method. We then introduce the real reservoir data set used to illustrate the method before explaining in more details the three different steps of the methodology: first, the reconstruction of full physics distances between models according to proxy distances; then the determination of the mathematical curve which best fits the oil production profiles; and last, the determination of the missing dynamic profiles by non-linear optimization. Finally, we present the results and discuss the ability to compute more accurate quantiles on production forecasting for uncertainty evaluation.

Problem definition

We consider a large set of reservoir models covering the range of possible reservoir geometries and heterogeneities .

(1)

Each of these models can be generated using alternative interpretations scenarios and Monte-Carlo sampling. This ensemble represents the known structural and parameter uncertainty about the reservoir. Reservoir production forecasts reflecting geological uncertainties should ideally consider all these possible models. However, to keep the computational cost realistic, multiphase flow simulation can only be performed on a subset of models deemed representative of the ensemble:

(2)

Finding these models is a problem per se, which calls for a priori ranking of all the models of based on their likely dynamic behavior. In a second step, the analysis of the production scenarios is usually based on production curve quantiles. The few selected models of are used to simulate the dynamic response and some specific quantiles (Q10, Q50 and Q90) are computed from the production curves at each time step (e.g., [START_REF] Subbey | Prediction under uncertainty in reservoir modeling[END_REF]).

The first limitation with this approach is that only the information from is taken into account and all the remaining models (from ) are discarded. The second limitation is that the quantiles are computed from very few data and therefore they are very sensitive to the model selection.

The goal of our methodology is to address these two challenges by: (1) using the results of flow simulation proxies on the models of ; (2) providing confidence intervals around the quantiles of the production curves, which reflect both the limited number of model realizations and the quality of the proxy.

Overview of the proposed methodology

Fig. 1 illustrates the overall methodology proposed in this paper. The general aim is to minimize the loss of information and the sensitivity of the quantile estimation. This is attained by estimating surrogate oil recovery profiles for all the models belonging to in order to account for the entire set of models. This allows to estimate the uncertainties more accurately while keeping the computational cost relatively low.

As illustrated on the top left of Fig. 1, we use a proxy for every model belonging to following the same philosophy as [START_REF] Josset | Local and Global Error Models to Improve Uncertainty Quantification[END_REF]. Each physically-based simulation proxy provides an approximate production curve , where refers to a given reservoir model and is a pseudo-time step. We use these data to compute a distance between each pair of models. The notion of distance between geomodels is presented by Suzuki and Caers (2008) to quantify the dissimilarity between the models. Here, we compute the squared distance between proxy curves as ,

where represents the proxy curve of the model at pseudo-time step , the proxy curve of the model at pseudo-time step and denotes the last time step of the simulation. Alternatively, proxies may consist of simple scalar measures , such as oil in place or pore volumes connected to producing wells. In these cases, the distance becomes:

,

The proxy distances are used to select representative models (the ensemble). For each of the selected models, the complete and accurate flow simulator is used to compute the dynamic flow response (red rectangle in Fig. 1), following the methodology of [START_REF] Scheidt | Representing Spatial Uncertainty Using Distances and Kernels[END_REF]. The distances between the accurate flow responses are computed as well.

Then the relations between the two sets of results are analyzed statistically. Using a simple analytical expression for the production profile (surrogate profile), the statistical analysis between the proxy and full physics response distances, we estimate the surrogate profiles for all the models belonging to (distance reconstruction and profile reconstruction in Fig. 1). This allows to generalize the workflow of [START_REF] Josset | Local and Global Error Models to Improve Uncertainty Quantification[END_REF] in the case of diverse types of proxies.

Finally, we estimate uncertainties on oil recovery by computing the production quantiles Q10, Q50 and Q90 using both the accurate responses obtained for the models belonging to and the surrogate profiles for the models belonging to (lower left corner in Fig. 1). All the steps are described in detail in the following sections.

Test case

To illustrate the proposed methodology and test both the workflow and its individual steps a real reservoir case study is used. The reservoir consists of turbiditic channel deposits and is bounded laterally by pelagic shales of negligible permeability. The reservoir layers have been deformed and affected by several faults (Fig. 2). The field was developed with four producing wells. Three water injection wells and an aquifer maintain the pressure during reservoir production. For confidentiality reasons, we cannot disclose reservoir location, dimensions and formation names.

In this case, the ensemble contains 163 corner-point reservoir models of cells describing the static reservoir uncertainties. This number of models corresponds to the subset of the initial 200 generated models for which the full physics flow simulations converged. We preferred to discard the models for which flow simulation did not converge in order to have a set of comparable models. Solving the convergence issues would require changing some numerical parameters and it could lead to differences between the model responses that would not be related to the geology. The grid geometry was generated from seismic interpretation and well tops. To reflect structural uncertainties, the reference structural interpretation was perturbed within bounds that reflect seismic time-to-depth conversion uncertainties and uncertainties in picking the reservoir layers on the migrated seismic image [START_REF] Abrahamsen | Bayesian kriging for seismic depth conversion of a multi-layer reservoir[END_REF]Thore et al., 2002). As a result, all grids in the set of realizations have distinct geometries. For the petrophysical modeling, we first simulated lithologies using a multiple point simulation technique, then we used sequential Gaussian simulation (SGS) within each facies for net-to-gross, effective porosity and co-located simulation for the permeability field. In the SGS, we used a proportion cube from a seismic attribute. For the simulations, we honored the target statistics inferred from the wells and spatial variability models (variograms) obtained from available data (wells) and analogs. Finally, the oil-water contact was sampled within a range of 20 meters and corresponding initial water saturations were simulated using again sequential Gaussian simulation (Fig. 2). The aim is to forecast field production for 20 years.

Proxy simulation was performed on all models. In another study [START_REF] Bardy | Proxy Comparison for Sorting Models and Assessing Uncertainty on Oil Recovery Profiles[END_REF], we investigated the quality of several proxies on the same ensemble of models. In this paper, we do not use an actual flow proxy, but simply consider the Stock Tank Original Oil In Place (STOOIP) as an indicator of reservoir production performance. STOOIP characterizes every model by an integer value using the following equation:

, ( 5
)
where: is the number of cells in the model, GRV (Gross Rock Volume) is the effective part of rock in the cell, NTG (Net To Gross) corresponds to the proportion of reservoir rock in the cell, is the water saturation in the cell and Bo the oil volume factor (volume ratio between reservoir and surface conditions). Although STOOIP is not a direct reservoir production proxy, it can be taken as such by our distance-based method. As expected theoretically and shown by [START_REF] Bardy | Proxy Comparison for Sorting Models and Assessing Uncertainty on Oil Recovery Profiles[END_REF], this proxy is less accurate than upscaling to provide a good approximation of model responses, but we chose it because it is very simple and quick to compute. As a result, it is often used in the oil industry to rank static reservoir models, even though the relationship between static and dynamic reservoir attributes can be complex. This choice was also made to evaluate to what extent our method can estimate production curves from a scalar static proxy. As this proxy only provides a single value per model, the distance was computed using Eq. ( 4).

The flow simulation forecast over 20 years were performed on all the stochastic models using the finite volume based software ECLIPSE® from Schlumberger®, configured with three fluid phases.

Note that another test case for the methodology is presented in [START_REF] Bardy | Intégration de modèles approchés pour mieux transmettre l'impact des incertitudes statiques sur les courbes de réponse des simulateurs d'écoulements[END_REF]. Injector (blue) and producer (green) wells trajectory are also visible.

Model selection

The selection of the subset must be as representative of as possible to keep the diversity into the responses and to allow an accurate uncertainty assessment on production. To achieve that, we use the Distance Kernel Method (DKM) proposed and described in detail by Scheidt andCaers (2007, 2009) as well as [START_REF] Josset | Local and Global Error Models to Improve Uncertainty Quantification[END_REF]. Here we provide only a brief summary of this technique. The interested reader is referred to the original papers for more details. The distance matrix between each pair of model can be computed by Eq. (3) or Eq. ( 4) depending on the retained proxy. Then the Multi-Dimensional Scaling algorithm (MDS) takes the proxy distance matrix as input and allows to plot the models in an abstract space of low dimensionality relatively to the number N of models. In this space, a set of representative models can finally be selected with a K-medoid algorithm (which is a clustering algorithm applied after a kernel transformation, see [START_REF] Kaufman | Clustering by means of medoids[END_REF]. This clustering algorithm finds models (medoids) which partition the model space so that the sum of distances between each model and the closest medoid is minimal. The main advantage of this method as compared to a classical K-means is that it is returning the medoids (models closest to the center of each cluster) which can be used directly, instead of the centers of each cluster which do not correspond to an existing model. The choice of is important because it should neither be too high (which will be too time consuming to simulate full physics) nor too small (which will not be sufficiently representative). The selected models belonging to on which flow simulations will be performed correspond to these medoids.

In our case study, we chose to use ten models ( ) as it is sufficient to correctly calibrate the method while minimizing the computational burden. As k-means, the k-medoid method solves a global optimization problem iteratively; as a result, different sets of models may be obtained

Missing Distance reconstruction

In the previous step, the distances between the proxy responses of the models was used to select the models belonging to the subset . The detailed and accurate flow simulation is performed only for the models belonging to , and results in an oil recovery profile with and . In most approaches using model selection, one would estimate the overall uncertainty on oil recovery by computing directly the quantiles from this limited set of profiles.

Here, we propose to go a step further and use all the information that was obtained by computing the proxies for all the models to improve the overall uncertainty quantification. The idea is that the distances between the proxy responses can be correlated with the distances between the real flow responses. Using that correlation and using in addition a simple expression for the shape of the profiles, it may be possible to reconstruct the variability of the unknown profiles reasonably well without solving the accurate flow problem and then use them to improve the uncertainty estimate. To do so, the first step consists in estimating the unknown distances between the true profiles from the distances between the proxy responses.

More precisely, from the accurate data, it is possible to compute the "true" squared distance ( ) between model responses:

.

For each pair ) of models belonging to , the distances between the proxy curves and the distances between the full physics curves are available.

Fig. 4a presents the cross-plot between full flow simulation distances and proxy distances obtained in this case; note that different selected models provide distinct distances. A linear regression between the two types of distances could be estimated but there is no theoretical reason to assume a linear relationship. A non-linear correlation is confirmed by our experiments (Fig. 4c), and therefore a piecewise linear correlation model on intervals is used instead. We choose these domains in order to contain approximately the same number of points ( or , Fig. 4a). We also take care of having enough points in each domain to be able to compute reasonable statistics.

As we have 10 simulated models, this means distinct differences which give 45 points, we choose 5 domains of 9 points each. For each domain, we use a simple linear model simulating the full flow simulation distance using a linear function of the proxy distance and a random error term :

(7)

The coefficients and of this regression are given by ordinary least squares.

Using equation ( 7), we simulate all the missing values of for all the pair of models belonging to . The random noise is randomly sampled from a normal distribution whose parameters are given by the linear regression model (centered Gaussian distribution) and corrected to ensure that only positive distances are obtained. Additionally, we force the coefficient to be 0 so that the regression passes through the origin in the first domain (a null distance corresponds to identical model response according to both proxy and full physics simulations). This piecewise regression allows for increasing variance with proxy distance, but may have discontinuities between domains and is not strictly consistent with regression theory.

In this work, we did not investigate more advanced regression models. Indeed, when applied to our test case (Fig. 4b), this method yields a set of acceptable possible full flow simulation distances which are consistent with the reference plot generated from the reference flow solution (Fig. 4c, [START_REF] Bardy | Proxy Comparison for Sorting Models and Assessing Uncertainty on Oil Recovery Profiles[END_REF] and which is sufficient for reasonable uncertainty estimates as shown below.

Note also that a perfect proxy, showing no error at all, would reproduce perfectly the true distances in the simulation space and would lead to a perfect straight diagonal line in Fig. 4a. However, a perfect proxy is very likely to be almost as computationally expensive as the accurate forward flow simulation and will have therefore no computational interest. Instead, a fast proxy will for sure display some discrepancies between the two distances as illustrated in Fig. 4a. The aim of the proposed method is therefore to model properly those errors for any kind of proxy and benefit from their numerical

Parametric model of oil production

Quantifying uncertainties about reservoir production calls for estimating oil production profiles for all the models belonging to (for which full physics simulations are unavailable). The profiles of total oil produced at field scale (FOPT) represent the cumulative oil production of the field at each reported time-step. In a standard workflow, these profiles are obtained from flow simulations, but proxies do not necessarily generate this type of curves. For example, when using STOOIP as a proxy, we can compute the distances between the models but not the FOPT curves.

This is why, in this work, we suggest estimating the FOPT profiles from the set of distances (Eqs. 6 and 7), the proxy simulations (Eq. 3) and the regression model (Eq. 7) using a parametric model (Eq. 8). The aim of the current section is to present the parametric model that we propose to use. In section 8, we will explain how the parameters of these curves are determined from the value of the distances.

A parametric model of the FOPT curve should respect three criteria imposed by the physics of oil recovery:

 It is a cumulative profile. The mathematical function must be monotonically increasingwith a positive or null derivative.  The first part of the profile corresponds to the maximum production rate in which liquid rates are constrained by the production installation (e.g., pipe diameter) so the function must be linear at the beginning and cross the origin.  The reservoir is considered at the beginning of its production, with no historical data and a constant number of production wells; after a while (time ), production declines so the function must have a strictly negative second derivative.

Rather than choosing between an exponential or hyperbolic analytical production decline model (Fetkovitch et al., 1996), we use a piecewise parametric function (Fig. 5a). It consists first in a linear part from the origin to a first control point corresponding to the beginning of the production decline, followed by a Cubic Hermite Spline from to a second final control point . As shown in Fig. 5a, the slope in Eq. ( 8) corresponds to the peak production rate and can be identified using the linear part on the simulated models during the early time steps. The final time corresponds to the duration of production so it is the same for all models. Therefore, five unknown parameters remain in Eq. ( 8) to fully define the analytical production profile :

--the time after which production starts declining. The cumulated production at this time is given by . --the total production at the final simulation time ( ).

--the norm of the tangent vector at the start of production decline ( ), We assume a gradual production decline so the direction of the tangent vector is the same as in the linear part, i.e., . In the Hermite polynomial, this norm controls how fast the production rate decreases after .

and -respectively the slope and the norm of the tangent vector at the end of simulation ( ), i.e.:

. The slope corresponds to the production rate at the final simulation time; and the norm controls how quickly this rate decreased in the late phases of production.

To test the relevance of this analytical production model, we considered the set of FOPT curves obtained by flow simulation on the reservoir models of , then we used least-squares minimization to fit the parameters of the analytical model (Eq. ( 8)). Results (Fig. 5b) suggest that the proposed function is able to adequately represent FOPT under our simulation settings (continuous field production without interruption). This step will also provide an ensemble of valid sets of parameters that we will use in the next part of the workflow.

Reconstruction of profiles using the proxy distance

Let us consider a reservoir model belonging to . The aim of this section is to reconstruct a possible FOPT curve for this reservoir model without running the flow simulation. We propose to do this by identifying the five parameters of the analytical function given in Eq (8) and described in the previous section. For this, we introduce a new squared distance between the analytical production profile for model and the known profile of each of the models of :

(9)

Ideally, the analytical profiles should be such that their distance to the known models of be equal to the distance simulated from the piecewise linear regression (Eq. 7). Note that we assume that is a reasonable approximation of the unknown distance between the true profiles.

Therefore, we propose to find the five unknown parameters of each profile (Eqs. 3-4) by minimizing the following objective function: (10) This minimization is a non-linear problem; it can be solved using the interior point method [START_REF] Nocedal | Numerical optimization[END_REF]. In practice, to help convergence, the curves are normalized so that the highest known profile is equal to 1.

Furthermore, to ensure an acceptable result, we prescribe bounds one some parameters and adde some non-linear constraints. To force the production decline to start before the end of the simulation period, we constrain to be in the interval . To reflect that production rate is lower at the end of the simulation than during peak reservoir performance, we constrain to be within . We also assume that the best model in terms of final production does not produce more than two times the production of the best simulated models:

. Last, and are bound within the interval to avoid sharp variations of the production curve between times and . In addition, we set a constraint to respect the physics of the simulation by forcing the second derivative of the mathematical model (Eq. ( 8)) to be negative (

).

In our study, we first applied this method to the set of ten reservoir models of , for which we know both the FOPT curve obtained by simulation and proxy results. For each model, we minimized Eq. ( 10) to estimate the analytical FOPT curve from the distances to the nine other models ( ). Fig. 6a shows that the ranking of the estimated analytical curves is the same as the reference ranking for most time steps (especially for long simulation times). However, some rank inversions are visible for intermediate time steps. Also, the linear part of the obtained production curves is often too short as compared to the true profiles. This could be due to the relatively small number of samples, or to the existence of local minima in the non-linear optimization. To reduce these errors we propose an optional kriging-based correction method described in Appendix.

We then ran the non-linear minimization of Eq. ( 10) for all the models of based on the distance to all ten models of set ( ). To speed up this optimization, we selected as starting point of the minimization the parameters of the known curve for which the objective function value was lowest.

Three individual comparisons between the estimated analytical profiles and the true profiles from the reference simulations are presented in Fig. 6b-d. In these cases, also, the analytical profiles are correlated to the reference ones, but display some variations. We think these deviations occur partly because of errors in the proxy-based distances and partly because distances are unable to uniquely characterize the shape of the production curves. Further research is needed to characterize these sources of error and to characterize the numerical properties of the minimization of Eq (10).

Production curve quantiles and confidence intervals

The above procedures produce a set of true production profiles (from flow simulation on models belonging to ) and estimated analytical profiles (from the minimization of Eq. ( 10) on all models of ). The corresponding curves are shown on Fig. 7a. For each time step, we create a Fig. 7b shows the results of quantile curves estimation for our case study. The curves estimated with our methodology are closer to the reference than the production profiles of the models in . This confirms experimentally that proxies integrated with the proposed method reduce uncertainty quantification errors as compared to uncertainties estimated only from a few representative models. However, the errors discussed in Section 8 are still visible when comparing the computed quantiles and the reference quantiles obtained from flow simulation on all models. In particular, the Q10 and Q50 curves underestimate the time when production starts declining, yielding more conservative production scenarios. To remove these systematic errors, we applied the correction described in Appendix. As shown on Fig. 7c, this correction introduces some non-admissible behaviors in individual curves. For instance, some cumulative production curves locally have a negative slope. However, these local inconsistencies are filtered by the quantile estimation process and the final quantile curve estimates (Fig. 7d) are closer to the reference curves.

In all cases, quantile curve estimations are prone to several sources of uncertainty. First, there is no guarantee about the representativeness of the subset of models selected for full flow simulation. This could ideally be addressed by increasing the number of models in until the output statistics converge. Because this can be computationally infeasible, another option is to repeat the k-medoid clustering process (Section 5) starting from different initial models. Indeed, the k-medoid algorithm is sensible to the choice initial models, as mentioned by [START_REF] Josset | Local and Global Error Models to Improve Uncertainty Quantification[END_REF]. This option is computationally expensive but can provide a practical way to assess confidence intervals around the quantile curves for a fixed number of models in .

A second source of uncertainty relates to the nature of the chosen proxy which, as the name suggests, only approximates the actual flow response. Indeed, more accurate proxies should in principle lead to better confidence in the quantile curve estimates. In our method, this proxy error is sampled during the missing distance reconstruction (Section 6). Resampling alternative distances from the piecewise regression model provides an efficient way to address this source of uncertainty without need for additional flow simulations.

In this study, we first applied our methodology 200 times. Due to the stochasticity of the k-medoid algorithm used in the DKM (Section 5), we were able to generate 200 different sets of quantile profiles (Q10, Q50 and Q90) and then computed the confidence interval around them with ten representative models (Fig. 8a). The confidence intervals obtained with our method do not always include the actual quantile curves generated from the reference simulations due to proxy error or workflow limitations. We suspect that the correction step described in the appendix could explain why these confidence intervals are too low. Indeed, this correction may artificially lower the profiles. Therefore, we chose to increase the number of representative models to twenty and ran the randomized workflow 200 times to obtain new confidence intervals (Fig. 8b). The results in this case are better in the sense that confidence intervals are slightly narrower and include the reference quantile curves. In contrast, we applied 200 random selections of representative models, and aggregated the true production curves to compute the associated confidence intervals (Fig. 8c). When our methodology is used, the confidence intervals are narrower due to the integration of proxy knowledge.

Conclusions

Estimating the uncertainties on the flow response of a reservoir by performing full physics flow simulation on a large ensemble of models implies often prohibitive simulation time. To overcome this problem, some approaches consider an approximation of the dynamic simulation response which is obtained faster by simplification of the physics of the simulator. The use of these proxies leads to a delicate issue of balancing acceleration with accuracy of the response. Another problem is that proxies are very often case-dependent; choosing the proxy which is best suited to a given case study can be a long procedure [START_REF] Bardy | Proxy Comparison for Sorting Models and Assessing Uncertainty on Oil Recovery Profiles[END_REF]. Other approaches consider selecting only a few models from the initial set and run the accurate flow simulation only on this small subset. Our approach aims at integrating the dynamic behavior of the whole set of models to capture the uncertainties, as also proposed by [START_REF] Josset | Local and Global Error Models to Improve Uncertainty Quantification[END_REF]. The method proposed in this paper allows for a larger class of proxies to be considered, but it calls for the explicit definition of analytical production profiles, which depend on the type of reservoir simulation output. From a reservoir management standpoint, it could be interesting to use other types of production curves such as water production or pressure, because it can also be important to adjust surface facilities. Looking at the well by well response would also be important as this is where proxies lack accuracy. Applying this methodology to any other kind of curve like water profile is possible but relies on the availability of compact analytical models for the considered reservoir response curve.

In this paper, the model selection was performed using the Distance Kernel Method [START_REF] Scheidt | Representing Spatial Uncertainty Using Distances and Kernels[END_REF] but other selection techniques could be considered. For field development studies, further research is also needed to assess how to address more complex development scenarios where multiple well locations are considered and well schedule changes though time.

Nonetheless, the results obtained show that the methodology can assess reservoir oil production uncertainty even with a very simple proxy. It would be interesting to assess the value of this method to accelerate history matching tasks in the same spirit as [START_REF] Scheidt | A multi-resolution workflow to generate high-resolution models constrained to dynamic data[END_REF] and Josset et al. (2015). It would also be interesting to test this methodology on a case study where the full flow simulation responses are more chaotic and with a less regular dispersion.

Fig. 1 -

 1 Fig. 1 -Proposed workflow. A few representative models are selected from proxy simulation (Section 5) and accurate simulations are performed. Using the correlation between proxy distances and accurate distances, all missing dynamic distances are determined (Section 6), then an analytic profile (Section 7) is fitted to determine the missing dynamic profiles (Section 8). Finally, oil recovery quantiles are computed (Section 9).

Fig. 2 -

 2 Fig. 2 -Our case study is made of several realizations (a-c: different geometry and facies distribution) in 3D view. Each realization is made of several petrophysical properties simulated successively (d: facies, e: porosity, f: permeability).

Fig. 3 -

 3 Fig. 3 -Examples of models selection. The red curves represent the 10 selected models while the green curves represent the other models. The vertical axis represents the cumulative oil production or Field Oil Production Total (FOPT)

Fig. 4 -

 4 Fig. 4-(a) Cross plot of full flow simulation distances vs. proxy distances for . Vertical lines represent domain borders. Red lines represent domain regression surrounded by the quantiles Q10 and Q90 of the residuals. (b) Cross plot of the reconstructed full flow simulation distances for (obtained by random draw as explained in the text) vs. proxy distances. (c) True cross plot of full flow simulation distance vs. proxy distance, used as benchmark

  Published in Journal of Petroleum Science and Engineering 176:116-125. https://doi.org/10.1016%2Fj.petrol.2019.01.035. This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License 11 where are the Hermite functions. Further symbols are defined below.

Fig. 5

 5 Fig. 5a) Graphical view of the analytical model and parameters of analytical oil recovery profiles. b) Result of leastsquares fit of the analytical model (red curves) to production profiles obtained by flow simulation (blue curves).

Fig. 6 -

 6 Fig. 6 -(a) Minimization of the analytical model parameters for the models belonging to , dashed curves correspond to the reference profiles obtained with a full physics flow simulator and the purple ones are obtained after minimization. (b-d) Examples of minimization procedure results. Red dashed curves correspond to the real response profile of models belonging to , green curves to the analytical curve of a model belonging to and blue curves to corresponding real profile (computed with full physics flow simulator and not known in operational context).

Fig. 7 -

 7 Fig. 7 -Final analytical production curves and corresponding quantiles (a) Analytical curves obtained after minimization applied to each model of (full green lines) and accurate curves obtained by flow simulation on models (dashed red lines). (b) Quantiles (Q10 in blue, Q50 in black and Q90 in red) computed with all models ( ) (reference, dashed lines); quantiles computed only with models belonging to (doted lines); quantiles computed with our methodology ( ) (full lines). (c) Minimized profiles after correction (see Appendix). (d) Quantiles updated with corrected profiles for our methodology (full lines).

Fig

  Fig. 8 -Quantiles dispersion over 200 iterations of the methodology. (a) With a subset of 10 models. (b) With a subset of 20 models. (c) With 10 models in but without using proxies from in the quantile computation.

Published in Journal of Petroleum Science and Engineering 176:116-125. https://doi.org/10.1016%2Fj.petrol.2019.01.035. This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License

Acknowledgments

The authors of this paper, would like to thank Total S.A. for funding this research and providing the case study, Celine Scheidt for providing the initial distance kernel method code, as well as Peter King (Imperial College) and Vicent Corpel (Total S.A) for their unconditioned help. We also thank Pedram Masoudi, and the other reviewers whose comments helped us to improve this paper.

Appendix: Error correction

We propose an optional correction stage to reduce the errors between the surrogate curves and the true curves. This correction makes the assumption that errors have zero mean, are smooth, and can be characterized only from the models of set . It is similar in spirit to the "gobal error model" introduced by Josset and Lunati ( 2013), but it considers production profiles instead of breakthrough curves and uses time-dependent weighting instead of a global distance-based weighting. This correction takes the form of a vertical shift map (Section 8), computed as follows:

-For every known model (belonging to ) the difference between the true profile and the reconstructed one is computed (Fig. 9a). -These differences are projected on the analytical curves in the plot space (Fig. 9b) -The differences are then interpolated for each time step using simple kriging with a 0 mean.

The juxtaposition of all the 1D corrections over all the time steps covers the entire domain FOPT versus time (Fig. 9c).

The final corrected production profiles are finally obtained by subtracting the error estimate from the reconstructed analytical profiles.