Toward an enhanced Bayesian estimation framework for multiphase flow soft-sensing

T o w ard an en h a n c ed Ba y esian estimati on framew ork for m ultiph a s e flo w soft-sensing Xiao dong Luo ∗ , Rolf J. Lorentzen, Andreas S. Storda l and Gei r Næ vdal In terna tiona l Researc h Instit ute o f Sta nv anger (IRIS) , 50 08 Berg en , Nor w ay Abstract In this w ork the authors study the multiphase flo w soft-sensing p roblem based on a previously established framework. There are three functional mo du les in this framew ork, namely , a transien t w ell flo w mo d el that describ es the r esp onse of certain physic al v ariables in a wel l, for instance, temp erature, v elocity and pressure, to the flo w rates en tering and lea ving the w ell zones; a Mark o v j ump pro cess t h at is desig n ed to captur e the p oten tial abru pt c han ges in the fl o w rates; and an estimation m etho d that is adopted to estimate th e underlying flo w rates b ased on the measurements from the physical sensors installed in the w ell. In the previous studies, the v ariances of the flo w rates in the Marko v jump pro cess are chose n manuall y . T o fill this gap, in the cur ren t w ork t w o automatic app roac hes are p rop osed in order to optimize the v ariance estimation. Through a numerical example, we sh o w that, wh en the estimation fr amew ork is used in conjun ction with these t wo pr op osed v aria n ce-estimat ion approac hes, it can ac hiev e r easonable p erfor- mance in terms of matc hing b oth the measur emen ts of the physica l sen sors and the true u nderlying fl o w rates. 1 In tro duction Smart wells are adv anced op eration facilities used in mo dern oil and ga s fields . T ypically , a smart we ll is equipped with do wnhole s ensors and inflow con trol v alv es (ICVs ). The ICVs are used to con trol the inflo w rates of the fluids in the w ell, whic h th us requires the information of the inflow rates. Nor ma lly , direct flo w ra te sensors are not oft en deplo y ed in a smart we ll, since they are exp ensiv e to use, y et the resulting measuremen ts ma y b e unreliable due to the high risk of sensor malfunction o r fa ilure. Therefore, in pra ctice it is more customary to install do wnhole sensors that collec t and transmit other types of dat a , e.g., pressure and temp erature. Based on the collected data, one can estimate the flow ∗ Corr esp ondi ng aut hor addr ess: Int er national Resear ch Institute of Stav a nger (IRIS), Thormøhlens Gate 55, 500 8 Bergen, Nor w ay . E -mail: xiao do ng.luo@iris.no 1 rates in the w ell. Suc h an exercise is often called “soft sensing” or “ soft metering” (see, for examples, Blo emen et al., 2006; de Kruif et al., 2008; Lesk ens et al., 2008; Loren tzen et al., 2010b; W rob el a nd Sc hiferli , 2009). Infor mat ion obtained f r om soft sensing can then b e used to supp o rt decision-making, e.g., c ho osing ICV op eration strategies fo r the purp o se o f pro duction optimization. In this sens e, soft sensing is essen tial for efficien t and profitable pro duction operat io ns. V arious metho ds hav e b een prop osed in the literature in order to tack le the soft sens- ing problem based on t he measuremen ts collected from down hole sensors. F or instance, Blo emen et al. (200 6) used the conv en tional extended Kalman filter (EKF) for soft sensing in gas-lift w ells; Gryzlo v et al. (2010); Lesk ens et al. (2008); Lo ren tzen et al. ( 2 010a,b) adopted the more recen tly emerged ensem ble Kalman filter (EnKF, see, for example, Aanonsen et al., 200 9; Ev ensen , 2 006; Nævdal et al., 2005) to estimate the flo w rates based on the temp erat ure and/or pressure measuremen ts. T o b etter a ddress the issue of non- linearit y in the pro cess of soft sensing, a more sophisticated metho d, called the auxiliary particle filter (APF, see, for example, Pitt and Sheph ar d, 19 9 9), is emplo y ed in a recen t w ork of Loren tzen et al. (2 014a). Apart from the aforemen tioned (approx imate) Ba y esian filters, the soft sensing pr o blem can also b e solved by alternative metho ds, fo r example, through a certain iterativ e algorithm in the in v erse problem theory (see, f o r example, Y oshiok a et al., 2 009). The curren t work adopts the same framew ork a s tha t dev elop ed in L o ren tzen et al. (2014a). There are three functional mo dules in this framew ork, including: (1) a transien t m ultiphase we ll flo w mo del, in contrast to the con ve ntional steady state mo dels used in man y other w orks (se e, f o r example, Azim, 2012). Here, the w ell flow mo del sim ulates the b eha viour of a s mart w ell with giv en flo w rates and other w ell conditions; (2) a Mark ov jump pro cess that describes how the flo w rates change with time. O ne reason for Loren tzen et al. (2014a) to ado pt this dynamical mo del with jumps is for its ability to capture the rapid c hanges in the flo w ra tes, whic h might tak e place due to the (sudden) changes of w ell conditions, for example, in case of water or gas breakthrough, or op erational c hanges of other we lls in the field, a nd so on (Loren tzen et al., 2014a); (3) an estimation metho d based on t he APF that estimates the underlying flow rates, based on av ailable measuremen ts from do wnhole sensors. Note that in general one can replace the APF by an y other suitable estimation metho d, e.g., those men tioned in the preceding paragraph. By using the aforemen tioned framew ork, satisfactory soft sensing r esults ar e obta ined in L o ren tzen et al. (20 14a). One o f the remaining issues in Lorentzen et al. (2014a) is that, when using the Mark ov jump pro cess, the v ariances of the pro cess noise are manually c hosen, whic h are th us not guarantee d to b e optima l in general situations. In the curren t w ork, we aim to fill this gap by prop osing tw o metho ds that are able to calculate the optimal v ariances based on certain optimalit y criteria. F or more inf o rmation a b out the transien t well flow mo del and t he APF, readers are referred to Lorentz en et al. (2014 a, 2010a,b) a nd Pitt and Shephard (1999), as w ell as the app endix in this w ork. This work is organized as fo llo ws. Section 2 outlines the problem that we are in terested in, and prop oses tw o metho ds to enhance the estimation framew ork established previously . Section 3 uses an example to illustrate the p erformance of the estimation framew ork that 2 is equipped with the prop osed metho ds. Section 4 concludes the work. F o r readers’ b etter comprehension of the framew ork used in this work, a short introduction to the APF is also pro vided in App endix A. 2 Metho dolog y 2.1 Problem statemen t Without lo ss of generality , suppo se that the (abstract) transien t w ell flow mo del is describ ed b y the follo wing mathematical equation y k = H k ( Q k ) + v k , (1) where y k is the nominal measureme nt con taining, for instance, the temp era t ur e and pres- sure at the time instant k ; H k is the w ell flow mo del; and Q k = [ Q k , 1 , · · · , Q k ,r ] T is the rates of inflo w or outflo w of fluid phases in r well zones. Since in practice the measure- men ts from the ph ysical sensors normally may con tain certain errors, an additional term, v k , is t h us in presence in order to accoun t for this p ossibilit y . In this w or k w e assume that v k follo ws the normal distribution N ( v k ; 0 , W k ), with zero mean and co v ariance W k . F or simplicit y , here w e further assume that W k is kno wn to us, o therwise it ma y a lso b e estimated by the metho ds presen ted below. The Marko v jump pro cess of the flo w r ate at the i -th w ell zone is g iv en b y Q k +1 ,i = θ k +1 ,i Q k ,i + u k ,i , i = 1 , 2 , · · · , r , (2) where θ k +1 ,i is a m ultiplier that follows a certain discrete pr o babilit y distribution, e.g., P r ( θ k +1 ,i = c j,i ) = P j,i ( j = 1 , 2 , · · · , N i ) for N i predefined constants c j,i (sp ecific to the w ell zone i ), and u k ,i the pro cess noise following the normal distribution N ( u k ,i ; 0 , σ k ,i ), with the v ariance σ k ,i to b e estimated. A few remarks are in order b efore w e pro ceed further. Firstly , in the con text of da t a assimilation nomenclature (see, for example, Ev ensen, 2006), Eq. (2) represen ts the dy- namical mo del whose states (the flo w rates) are to b e estimated, while Eq. (1) is the asso ciated observ ation op erator that describ es the relation b etw ee n the mo del states and the observ a tions. Ba sed on this p oint of view, the noise term u k ,i in Eq. (2) can b e considered a s one of the sources of mo del errors. In general, the distribution of u k ,i is unkno wn. Here w e assume that u k ,i follo ws a Gaussian distribution with zero mean but unkno wn v ariance σ k ,i . Later we will discuss ho w to optimize the c hoice of σ k ,i (under a certain optima lity criterion) in o rder to improv e t he p erfo rmance of the Bay esian da ta assimilation algorithm in use (see, for example, Appendix A). Secondly , with the discrete random v ariable θ k +1 ,i and the Gaussian noise u k ,i in Eq. (2 ), it can b e ve rified that the corresp onding flow rates Q k +1 ,i is described b y a Gaussian mixture mo del (GMM). In data assimilation practices, using a GMM to appro ximate the probability density function (p df ) of a r a ndom v a r iable has certain adv an tages. F or instance, from a theoretical p oint of view, 3 the GMM is one t yp e of the k ernels in k ernel densit y estimation (KD E) theory a nd can b e used to a pproximate a generic non- Gaussian pdf (Silverman, 1986) 1 . In addition, using a Gaussian noise u k ,i in Eq. (2) is also con ve nient for our t heoretical dev elopmen t later. F or instance, this leads to a con v ex optimizatio n problem in Section 2 .3 .2 that can b e more con v enien tly solved, in con trast to p ot entially non-conv ex for ms with other distribution c hoices for u k ,i . In Loren tzen et al. (2014 a) the authors adopted an alternativ e no nlinear and non- Gaussian data assimilation algorithm, namely , the auxiliary part icle filter (APF, see Pitt and Shephard, 1999), to estimate t he flo w rates ba sed on av ailable measuremen ts. The difference betw een the APF and the con ve ntional sequen tial imp ortance re-sampling particle filter (SIR-PF, see, for example, G o rdon et al., 1993) is mainly at the re-sampling step. In t he SIR-PF, the re-sampling is conducted based on the most recen t p osterior w eigh ts of the particles, while in the APF, measuremen ts ahead (in time) of the particles to b e re-sampled are also tak en in to accoun t, suc h that after re-sampling the generated particles ma y tend to matc h the “future” measuremen ts b etter. F or brevit y , w e outline b elow the main pro cedures at the re-sampling step of the APF. Readers are r eferred to the app endix in t his w ork and Pitt and Shephard (1999) for more details. In the contex t of flo w rate estim at io n, the r e- sampling s tep in the APF is as fo l- lo ws. Supp ose that Q j k = ( Q j k , 1 , Q j k , 2 , · · · , Q j k ,r ) T is t he j -th sample in the set { Q j k } n j =1 , and w j k is the posterior w eigh t asso ciated with Q j k . F or eac h Q j k , a sample θ j k +1 ≡ ( θ j k +1 , 1 , θ j k +1 , 2 , · · · , θ j k +1 ,r ) T is also dra wn base d on the dis tributio ns of P r ( θ k +1 ,i ) ( i = 1 , 2 , · · · , r ). Let Ω j k +1 = ( θ j k +1 , 1 Q j k , 1 , θ j k +1 , 2 Q j k , 2 , · · · , θ j k +1 ,n Q j k ,n ) T , then one can compute an in termediate weigh t µ j k +1 ∝ w j k N ( y k +1 ; H k +1 ( Ω j k +1 ) , W k +1 ), where N ( y k +1 ; H k +1 ( Ω j k +1 ) , W k +1 ) denotes the p df of the Gaussian distribution of y k +1 with mean H k +1 ( Ω j k +1 ) and cov ari- ance W k +1 . Based on µ j k +1 , one applies SIR to the indices j = 1 , 2 , ...n and obta ins a new set of indices j s , s = 1 , 2 , · · · , n to pic k up samples from the set { Ω j k +1 } n j =1 . The in- dexed samples { Ω j s k +1 } n j =1 has equal we ights and is used t o g enerate new flow rates samples Q s k +1 ,i through Q s k +1 ,i = Ω j s k +1 ,i + u k +1 ,i ( i = 1 , 2 , · · · , r a nd s = 1 , 2 , · · · , n ). Finally , the w eigh ts w s k +1 asso ciated with t he sample Q s k +1 = ( Q s k +1 , 1 , Q s k +1 , 2 , · · · , Q s k +1 ,r ) T is computed b y w s k +1 ∝ N ( y k +1 ; H k +1 ( Q s k +1 ) , W k +1 ) / N ( y k +1 ; H k +1 ( Ω j s k +1 ) , W k +1 ). In the course of generating samples Q s k +1 through the form ula Q s k +1 ,i = Ω j s k +1 ,i + u k +1 ,i , the v ariances σ k +1 ,i of the Gaussian random v ariables u k +1 ,i are c hosen man ually in Lorentze n et a l. (2014a), whic h implies that in general there needs to b e a trial- and- error pro cess, and that the optimality of the final c hosen v alues ma y not b e clear. In this w ork we consider t wo automatic approa c hes that can b e used to optimize, in the sense to b e sp ecified lat er, the c hoices of σ k +1 ,i based on a v ailable measuremen ts. In the first approac h, the optimal v ariances are considered constant in time, and their estimation is obtained by minimizing a cost f unction with resp ect to the collected measuremen ts in a 1 In accor dance with this viewp oint, o ne may c ho ose other t yp es of kernels in appr oximation, e.g., by letting u k,i in E q. (2) follow s ome alter na tive distr ibutions. An inv estigatio n in this asp ect will, howev er, not b e car ried out in the current work. 4 giv en time in terv al. This ty p e of estimation approac h can b e classified as a “fixed in terv al smo other” (Simon, 2006), therefore fo r reference it is called the “fixed in terv al smo othing approac h” in the curren t w ork. In the second approac h, the optimal v aria nces are up dated with time, while the estimation is carried out b y minimizing a cost function with respect to the measureme nts up to (and including) certain fixed time steps ahead of the v ariances to b e estimated. This cor r espo nds to a “fixed lag smo other” (Simon, 2006), hence is called the “fixed lag smo othing a pproac h”. 2.2 The fixed int erv al smo othing approac h The fixed in terv al smo o thing approa c h is based on the idea in L o ren tzen et al. (2014 a). F or nota t io nal conv enie nce, let Y o 0: K ≡ { y o 0 , y o 1 , · · · , y o K } denote the collection of a set of observ ed measuremen ts y o • from time instant 0 to K , with K b eing fixed, and σ ≡ ( σ 1 , · · · , σ r ) T a vec tor con taining the v ariances of the Mark ov jump pro cess in differen t w ell zones. Our ob jectiv e is to estimate a constan t v ariance (v ector) σ that minimizes a certain cost function, as will b e explained b elow. T o construct the cost function, w e first denote b y p ( Y o 0: K | σ ) the j o in t proba bilit y densit y function (p df ) of Y o 0: K conditioned on σ . Then the optimal σ opt is tak en as the one that solv es t he follo wing minimization problem, i.e., σ opt = min σ ≥ 0 C ( σ ) , C ( σ ) ≡ − log p ( Y o 0: K | σ ) = − K X k =1 log p ( y o k | Y o 0: k − 1 ; σ ) , (3) where in the second line of Eq. ( 3 ) the identit y p ( Y o 0: K | σ ) = K Q k =1 p ( y o k | Y o 0: k − 1 ; σ ) is used. Note that, in Eq. (3), C ( σ ) is the sum of the terms − log p ( y o k | Y o 0: k − 1 ; σ ) at individual time instan ts, in whic h w e ha v e a ssumed that σ is constan t ov er time 2 . Therefore, in order to optimize σ , the observ ations in a fixed time in terv al (e.g., from 0 to K in Eq. (3)) need to b e used. In the litera t ur e, the corresp onding da ta assimilation sc heme is often called fixed in terv al smo other (see, for example, Simon , 2006), in whic h one needs to w ait un til all the observ ations in the sp ecified time in terv al arriv e b efore starting t he estimation. Because of the relatively long w aiting time (e.g., 100 min utes for our case study in Section 3), the fixe d in terv al smo ot hing approac h is often used for offline analysis . F urther note that p ( y o k | Y o 0: k − 1 ; σ ) = Z p ( y o k | Q k ) p ( Q k | Y o 0: k − 1 ; σ ) d Q k , (4) where p ( y o k | Q k ) is the conditiona l p df o f y o k on a giv en flow rates v ector Q k , and can b e determined through Eq. (1); p ( Q k | Y o 0: k − 1 ; σ ) is the prior (or predictiv e) p df of the flo w 2 In principle, o ne may choo se to let σ change with time (e.g., by replacing p ( y o k | Y o 0: k − 1 ; σ ) with p ( y o k | Y o 0: k − 1 ; σ k )) in E q. (3). In doing s o, how ever, the search spa ce of the minimization pro blem is larger a nd the optimization pr o cess b eco mes more exp ensive. 5 rates Q k , conditioned o n σ and the past measuremen ts up to (and including) time instan t k − 1. F or a give n σ , p ( Q k | Y o 0: k − 1 ; σ ) can b e obtained from a Bay esian filter. F or instance, in the con text o f particle filtering, p ( Q k | Y o 0: k − 1 ; σ ) ≈ n X h =1 w pr h δ ( Q k − Q k ,h ) , (5) where { Q k ,h } n h =1 is a set of n sample s (often called “particles”) of t he flow rates Q k ; { w pr h } n h =1 is a set of a priori weigh ts asso ciat ed with the particles, b efore the measuremen t y o k is incorp orated in estimation; and δ ( • ) represen ts the Dirac delta function, whose v alue is + ∞ at the origin, and 0 elsewhere. Substituting Eqs. (4) and (5) in to Eq. (3) , one obtains an approximate cost function, in terms of ˜ C ( σ ) ≈ − K X k =1 log n X h =1 w pr h p ( y o k | Q k ,h ) ! . ( 6 ) The impact of σ on the cost function is implicitly made through Eq. (2), whic h influences the generated particles Q k ,h . With the approximate cost function, a suitable optimization algorithm can then b e employ ed to optimize the v alue of σ . 2.3 The fixed lag smo othing approac h In the fixed lag smo othing approac h, the v ariance σ k c hanges with time, instead of b eing constan t. The amount of measuremen ts used to estimate σ k is , in most cases, also different from that in the fixed in terv al smo othing approach. In general, if one w ants to estimate σ k , it is the collection of the measuremen ts Y o 0: k + l ≡ { y o 0 , y o 1 , · · · , y o k + l } ( k + l ≤ K ) that is used, with l b eing a fixed (often small) p ositive inte ger (but note that the sum k + l c hanges with k ). Compared with the fixed interv al smo othing approach, its fixed la g coun terpart tends to b e more suitable for real time, online estimation ta sks, since it only needs to w ait for the l - step-ahead observ ations (i.e., observ ations from time instan t k + 1 to k + l ) to arriv e b efore the estimation of σ k starts, rather than all the o bserv ations fro m time instant 0 t o K a s in the fixed in terv al smo othing appro a c h. In the curren t study , l = 1 is c hosen in the analysis b elow, so t hat the metho d is also called lag -1 smoo t hing approac h hereafter 3 . In general, one ma y also consider long er lags, whic h, ho w eve r, is b ey ond the scop e of the curren t w ork. The optimization of σ k is done in a w ay similar to the exp ectation- maximization (EM) algorithm (D empster et al., 1 977). The EM algorithm is an iterativ e metho d and consists of t wo steps at eac h iteration. A t the exp ectation step (E-step), o ne constructs a cost function that can b e interpreted as the expectatio n of a certain quantit y . The cost function dep ends on the parameter(s), sa y σ j k , that will b e estimated at the iteratio n step j , a nd 3 Under this ch o ice, the w aiting time for the one-step-a head obser v ation is 2 minutes for our case study in Section 3, which appea rs more practica l, in co n tra st to 100 min utes in the fixed interv al smoothing approach. 6 Figure 1: An o ut line of the workflo w in the EM algo r ithm. it is also conditioned on the know n para meter v alue, sa y σ j − 1 k , that is obtained at the previous ( j − 1)-th iteration step. Give n t he cost function, the maximization step ( M-step) then tak es place in order to estimate an optimal v alue σ j k that maximizes the cost f unction v alue. Fig. 1 outlines the workflo w o f t he EM a lgorithm. A t the initialization step one sp ecifies the initial parameter v alue, sa y σ 0 k . One then constructs an exp ectation-t yp e cost function tha t is dep enden t o n σ 1 k and conditioned on σ 0 k . By maximizing the cost function, one obtains the optimal parameter v alue of σ 1 k . One then c hec ks if t he stopping condition is met. If not, a new iteration starts, in whic h one constructs a new cost function that is dep enden t on the parameter σ 2 k to b e estimated, and conditioned on the optimal v alue of σ 1 k obtained at the previous step, and so o n. The stopping condition in our implemen tatio n is either (a ) the relativ e Euclidean nor m c hange of tw o consecutiv ely estimated parameters, in term of k σ j k − σ j − 1 k k 2 / k σ j − 1 k k 2 , is lo w er than a pre-chos en threshold; or (b) the maximum iteratio n n um b er is reache d. T he estimate of σ k is then t a k en as the final one of the iterat io n pro cess. In what follows w e describe the E-step and M-step with more details. 2.3.1 The exp ectation step Without loss o f generalit y , supp ose that the v ariance σ k at time instant k needs to b e estimated. Also supp ose that an estimate σ j − 1 k is already obtained at the ( j − 1 )-th iteration of the EM algorit hm, so that the E-step aims to construct a cost function that dep ends on σ j k . In line with the c hoice l = 1 , the measuremen ts Y o 0: k +1 ≡ { y o 0 , y o 1 , · · · , y o k +1 } are inv olv ed in constructing the cost function, whic h is giv en b y E ( σ j k | σ j − 1 k ) ≡ Z log f  y , Q k , Q k +1 , θ k +1 | σ j k , Y o 0: k  f  y , Q k , Q k +1 , θ k +1 | σ j − 1 k , Y o 0: k  × δ ( y − y o k +1 ) d y d Q k d Q k +1 d θ k +1 . (7) In Eq. (7), y is a dummy v ariable with resp ect to the measuremen t y o k +1 . Since y o k +1 is already obta ined, the Dirac delta function δ ( y − y o k +1 ) is thus intro duced; θ k +1 ≡ 7 ( θ k +1 , 1 , · · · , θ k +1 ,r ) T is the a ugmen ted v ector of the m ultipliers θ k +1 , 1 , · · · , θ k +1 ,r in the Mark o v jump pro cess Eq. ( 2 ). Eq. (7) can b e interpre ted a s the exp ectation of the log lik eliho o d of y o k +1 , Q k , Q k +1 and θ k +1 (conditioned on σ j k and Y o 0: k ), with resp ect to t he join t p df of y o k +1 , Q k , Q k +1 and θ k +1 (conditioned on σ j − 1 k and Y o 0: k ). Let σ • k represen t either σ j − 1 k or σ j k , then one can fa ctorize the join t p df f  y o k +1 , Q k , Q k +1 , θ k +1 | σ • k , Y o 0: k  as follows: f  y o k +1 , Q k , Q k +1 , θ k +1 | σ • k , Y o 0: k  = f  y o k +1 | Q k +1  f ( Q k +1 | Q k , θ k +1 , σ • k ) f ( Q k | Y o 0: k ) f ( θ k +1 ) . (8) In Eq. (8), the conditional p df f  y o k +1 | Q k +1  can b e calculated based o n the t r ansien t well flo w mo del Eq. (1); the conditional p df f ( Q k +1 | Q k , θ k +1 , σ • k ) and the p df (with discrete v ariables) f ( θ k +1 ) a r e determined in the Marko v jump pro cess Eq. (2); the conditional p df f ( Q k | Y o 0: k ) is the p osterior p df of the flow rates, conditioned on the past measuremen t Y o 0: k up to (and including) time instan t k . Similar to the situatio n in the fixed in terv al smo othing approac h (see Eq. (5)), in the context of particle filtering, f ( Q k | Y o 0: k ) can also b e approximated by the empirical posterior p df of the samples of t he flo w rates. Based on the factor izat io n in Eq. (8), the log lik eliho o d term in Eq. (7) can b e re- written as log f  y o k +1 , Q k , Q k +1 , θ k +1 | σ j k , Y o 0: k  = lo g f  Q k +1 | Q k , θ k +1 , σ j k  + log  f  y o k +1 | Q k +1  f ( θ k +1 ) f ( Q k | Y o 0: k )  . (9) The second term on the right hand side of Eq. (9) is indep enden t of σ j k . Therefore, at the M-step this term can b e ignored for the purp ose of optimization. 2.3.2 The maximization st ep Com bining Eqs . (7) - (9 ), at the maximization step one needs to solve the f ollo wing optimization pro blem arg max σ j k ˜ E ( σ j k | σ j − 1 k ) ≡ Z log f  Q k +1 | Q k , θ k +1 , σ j k  f  y o k +1 | Q k +1  f  Q k +1 | Q k , θ k +1 , σ j − 1 k  × f ( θ k +1 ) f ( Q k | Y o 0: k ) d Q k d Q k +1 d θ k +1 , (10) with the simplified cost func tio n ˜ E ( σ j k | σ j − 1 k ) by discarding the irrelev an t term in Eq. (9). Since in g eneral it is difficult to obtain the analytical for m of the integral in (10), it is customary in practice to appro ximate that integral t hr o ugh certain Mon te Carlo appro ximations. Sp ecifically , let { Q k ,h } n h =1 b e a set of flo w rates samples at time instan t k , and { w post h } n h =1 the asso ciated po sterior we ights, then f ( Q k | Y o 0: k ) ≈ n X h =1 w post h δ ( Q k − Q k ,h ) . (11) 8 Similarly , let { θ k ,l } m l =1 b e a set of m ultiplier samples, then f ( θ k +1 ) d θ k +1 ≈ 1 m m X l =1 ˜ δ ( θ k +1 − θ k +1 ,l ) , (12) where ˜ δ is the Kroneck er delta function whose v alue is 1 at the origin a nd 0 elsewhere. The constan t factor 1 /m is not essen tial in optimization and is t hus dropp ed lat er. Giv en a flo w rate sample Q k ,h ≡ ( Q k ,h, 1 , · · · , Q k ,h,r ) T and a multiplier sample θ k +1 ,l ≡ ( θ k +1 ,l, 1 , · · · , θ k +1 ,l,r ) T , f  Q k +1 | Q k ,h , θ k +1 ,l , σ j − 1 k  follo ws a normal distribution with the mean Q k ,h · θ k +1 ,l ≡ ( θ k +1 ,h, 1 Q k ,l, 1 , · · · , θ k +1 ,h,r Q k ,l,r ) T , and the v ariances b eing σ j − 1 k ,i for i = 1 , · · · , r . F or Monte Carlo approximation, note that f  Q k +1 | Q k ,h , θ k +1 ,l , σ j − 1 k  d Q k +1 = f  Q k +1 | Q k ,h , θ k +1 ,l , σ j − 1 k  N ( Q k +1 ; µ , Σ ) N ( Q k +1 ; µ , Σ ) d Q k +1 , where N ( Q k +1 ; µ , Σ ) is the multiv ariate normal distribution with the mean b eing µ and the co v ariance b eing Σ . With these, one draw s a set of s samples, say { Q k +1 ,q } s q =1 , from N ( Q k +1 ; µ , Σ ), suc h that f  Q k +1 | Q k ,h , θ k +1 ,l , σ j − 1 k  ≈ s X q =1 δ ( Q k +1 − Q k +1 ,q ) f  Q k +1 | Q k ,h , θ k +1 ,l , σ j − 1 k  N ( Q k +1 ; µ , Σ ) . (13) It is natural to let µ b e the sample mean o f Q k ,h · θ k +1 ,l (with differen t com binations of Q k ,h and θ k +1 ,l ); o n the other hand, in our implemen tation w e let t he v ariances of Σ b e 3 times larger than t he empirical v ariances o f Q k ,h · θ k +1 ,l , suc h that it is (mor e) lik ely to ha v e certain samples a mong { Q k +1 ,q } s q =1 that are relative ly fa r a w ay fro m the sample mean of Q k ,h · θ k +1 ,l , ev en in the case that the empirical v aria nces of Q k ,h · θ k +1 ,l themselv es a re relativ ely small. Substituting Eqs. (11) - (13 ) in to (10), one obtains a Monte Carlo approximation of ˜ E ( σ j k | σ j − 1 k ), in terms of ˆ E ( σ j k | σ j − 1 k ) = X h,l, q Ω j h,l, q log f  Q k +1 ,q | Q k ,h , θ k +1 ,l , σ j k  , Ω j h,l, q = w post h f  y o k +1 | Q k +1 ,q  f  Q k +1 ,q | Q k ,h , θ k +1 ,l , σ j − 1 k  N ( Q k +1 ,q ; µ , Σ ) . (14) In order for σ j k to solve the maximization problem, a nece ssary condition is that ∂ ˆ E ( σ j k | σ j − 1 k ) ∂ σ j k = 0 . (15) Solving (15 ), one o btains the follo wing iterat io n form ula σ j k ,i = P h,l, q Ω j h,l, q ( Q k +1 ,q ,i − θ k +1 ,h,i Q k ,l,i ) 2 P h,l, q Ω j h,l, q , for i = 1 , 2 , · · · , r , (16) 9 whic h is an (iterativ e) w eigh ted squared error estimation in ligh t of the Mark ov jump pro cess Eq. (2). Giv en Eqs. (14) and (16), the iterat io n is done in the follo wing w ay: the v ariance estimate σ j − 1 k obtained a t the ( j − 1)-th iteration is used to calculate the w eigh ts Ω j h,l, q according to Eq. (14); then based on Eq . (16), the v ariance estimate is up dated to σ j k = ( σ j k , 1 , · · · , σ j k ,r ) T , which can then b e used to calculate the w eigh ts Ω j +1 h,l, q at the next iteration, and so on. Note t ha t b y Eq. (1 6), the v ariance estimates are p ositiv e definite (almost surely). In addition, when calculating the (itera t ive) we ights in Eq. (16), t he comp onen t w post h f  y o k +1 | Q k +1 ,q  / N ( Q k +1 ,q ; µ , Σ ) only needs to b e calculated once for all, while it is only the comp onen t f  Q k +1 ,q | Q k ,h , θ k +1 ,l , σ j − 1 k  that needs to b e iterat ively up dated. This feature significan tly reduces t he computational cost of the iteration pro cess. Finally w e note that, since f  Q k +1 | Q k ,h , θ k +1 ,l , σ j k  follo ws a normal distribution, it is easy to verify that the Hessian ∂ 2 ˆ E ( σ j k | σ j − 1 k ) /∂ 2 σ j k of the approx imate cost function ˆ E ( σ j k | σ j − 1 k ) is negativ e definite, which implies that the iteratio n formula Eq. ( 1 6) w ould in- deed lead to a v ariance estimation that maximizes the approx imate cost function ˆ E ( σ j k | σ j − 1 k ). 3 Example In this section w e illustrate, through an example, the p erformance of the APF with the optimized v ariance estimates obtained from the fix interv al/lag s mo othing a pproac hes. This example is taken fr o m the case study 2 in Lorentz en et al. (2014a). Note that, due to some recen t technic al up dat es, the transien t w ell flow mo del we are curren tly using ma y b eha v e sligh tly differen t than the one prev iously presen ted in Loren tzen et al. (201 4 a). In what fo llo ws w e consider tw o scenarios. In the first one, the w ell flo w mo del is run with a resolution o f 50m (in measured depth, md for short) in order to pro duce the syn thetic observ at io ns (pressure and temp erature), while in the course of data assimilation later, the we ll flow mo del is run with the same resolution (50m). F or distinction lat er, we call this p erf e ct observation sc enario (POS). In the second scenario, the w ell flo w mo del is run with a finer resolution of 5m t o pro duce t he syn thetic observ ations, while in the course of dat a assimilation, the well flow mo del is run with a coarser resolution of 50m. This in tro duces some additional observ ation errors to the observ ation system, apar t fro m the noise term v k in Eq. (1). F or this reason, w e call this imp erfe ct observation sc enario (IOS). In addition, in a few cases later w e also consider the p ossibilit y that the distribution of the observ ation error v k in data assimilation may mismatc h the o ne used to pro duce the syn thetic observ atio n. Note that our purp ose here is to examine the p erforma nce of the APF in the IOS, a nd no bias correction pro cedure is in tro duced in data assimilation to accoun t for the imp erfection of the observ at ion system. All the exp erimen ts b elo w are rep eated fo r a few times ( each time with a differen t seed in the random num b er generator), and similar r esults are o bserv ed. F or conciseness, in what fo llows we o nly presen t the results with respect to one o f the ra ndom seeds. 10 Figure 2 : W ell configur a tion. 3.1 The p erfect observ at ion scenario (POS) In the POS, the geometry o f the w ell is the same as that in Loren tzen et al. (20 14a). As sho wn in Fig. 2, the we ll has a v ertical part and a horizon tal part, connected through an inclined part in-b etw ee n. The inclined part has a curv ature of 0.11 degrees/m, starting from the v ertical depth of 1000 m and ending up with the vertical depth of 1500 m. In the horizon tal part, there are tw o influx zones: one lo cates b et w een md 3500 - 355 0 m (lab elled a s Z1), and the other betw een 3950 - 4000 m (lab elled as Z2 ) . In Z1 there is only gas influx, while in Z 2 there is only oil influx. There are tw o sensors (gauges) installed in the we ll, with the md b eing 3475 m and 392 5 m, resp ectiv ely (b oth gauges are placed 25 m ab o v e the influx zone). Thes e sensors collect the temp era t ur e and pressure at their lo cations. In our experimen tal settings, the reserv oir temp erature and pressure in Z2 a re 335.5 K and 1 . 5 × 10 7 P a, resp ectiv ely , while in Z1, they are 3 25.5 K and 1 . 4 × 10 7 P a, resp ectiv ely . F or mor e details, readers are referred to Loren tzen et al. (201 4a). Fig. 3(a) depicts the true flo w ra tes of the ga s and o il influx in Z1 and Z2, resp ectiv ely . After a spin-up p erio d, the w ell flow mo del starts t he sim ulation from time 10000 s until time 1 6 000s. During this p erio d the gas flow rate undergo es tw o abrupt c hanges: one o ccurs at time 106 00s, when the flow ra te increases from 2 kg/ s to 3 kg/s; t he other tak es place a t time 1280 0s, when the flo w r ate reduce s from 3 kg/s to 1 kg/ s. The oil flow rate also exp eriences one jump fro m 10 kg/ s to 15 kg/s, whic h happ ens a t time 10600s. In accordance with the ab o v e settings, Fig. 3(b) also show s the pressure and temp eratur e recorded by the gauges pla ced near Z 1 and Z2, resp ectiv ely (s ee Fig. 2). These data are collected ev ery 120s, and are con taminated with the Gaussian noise N ( v k ; 0 , W k ) (cf. Eq. (1)). Here W k is a diag o nal matrix, meaning that we assume the observ a tion errors of pressure and temp erature are uncorrelated with each other. In the exp erimen t, the standard deviations of the measuremen t noise are 0 . 02% of the magnitudes of the measured pressure and temperature data. 11 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 0.5 1 1.5 2 2.5 3 3.5 True profile of the gas flow rate Time [s] Zonal flow rate [kg/s] 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 9 10 11 12 13 14 15 16 17 18 True profile of the oil flow rate Time [s] Zonal flow rate [kg/s] (a) T rue flow rates 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3475 m Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 320 330 340 350 360 Sensor @ 3475 m Temperature [K] Time [s] 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3925 m Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 338 340 342 344 346 348 Sensor @ 3925 m Temperature [K] Time [s] (b) Measured pressure a nd temp era tur e Figure 3: (a) T rue gas and oil flow rates in Z1 and Z2 in the POS; (b) Measured pressure and temp eratur e from the gauges in the POS. The configuration of t he Mark ov jump pro cess is the same as that in Loren tzen et al. (2014a). Sp ecifically , for b oth t he gas and oil flo w rates, the m ultipliers θ in Eq. (2) tak e v alues from the same finite set { 0 . 5 , 0 . 75 , 1 , 1 . 25 , 1 . 5 } , while the corresp onding probabilities of taking those v alues are sp ecified b y the set { 0 . 1 , 0 . 1 , 0 . 6 , 0 . 1 , 0 . 1 } . The APF is initialized with the true flo w rates at time 10000s m ultiplied b y some samples o f θ dra wn at random. With the ab ov e said, in what follo ws w e presen t the n umerical results of the APF, with the v ariances of the Mark o v jump pro cess b eing manually t uned as in Lo ren tzen et al. (2014a), and those o bt a ined from the fix interv al/lag smo othing approaches . In all these cases, the sample size o f the pa r t icles is 500. 3.1.1 Results with the man ually tuned v ariances F or reference, w e first rep ort the results o f the APF with manually tuned v aria nces. F ol- lo wing Lor en tzen et al. (2 0 14a), we let the v ariances of the gas and oil flow rates b e b oth at 0.5. These man ually c hosen v ariances are used in the Marko v jump pro cess Eq. (2), based on whic h the flow rate samples (part icles) are draw n at random. The generated samples are then take n as the inputs to the w ell flow mo del, and the corresp onding outputs of temp erature and pressure a re compared to those recorded b y the gauges. Based on the data misfit b et we en the simulated temp erature and pressure and those from the gauges, the APF will adjust the relativ e w eigh ts of the particles, so that t he w eigh ted mean of the particles ma y matc h the observ ed temperat ure and pressure b etter. Since a set of samples is used, it also allows one to ev aluate the uncertain ties in the pro cess of es timatio n. Fig. 4 shows the estimated flow rates (pa nel (a)) and the corresp onding sim ulated temp erature and pressure (panel (b)). As one can see in Fig. 5(b), the sim ulated mea- suremen ts matc h those recorded b y the gauges well. On t he other ha nd, b efor e the second abrupt c hange of the gas flow rate happ ens (at time 12800s), the estimated gas and oil flow 12 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 0.5 1 1.5 2 2.5 3 3.5 Estimated parameter p t , p a ± 2 σ Time [s] Q 1 [kg/s] 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 9 10 11 12 13 14 15 16 17 18 Estimated parameter p t , p a ± 2 σ Time [s] Q 2 [kg/s] (a) Flow rates 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3475 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 320 330 340 350 360 Sensor @ 3475 m, d o , d a Temperature [K] Time [s] 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3925 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 338 340 342 344 346 348 Sensor @ 3925 m, d o , d a Temperature [K] Time [s] (b) Pr e s sure and temp era ture Figure 4: P erforma nce o f the APF with ma nually tuned v ariances. (a) Profiles of the gas and oil flo w rates in Z1 and Z2. The true ones are in red, and the estimated ones are in blue (solid). Also plotted in blue (dash-dotted) are the flow ra tes that are ± 2 standa r d deviations aw a y fro m the estimated ones; (b) Pressure and temp erature. Those recorded b y the gauges are in red (dotted), while those sim ulated o nes based on the estimated flow rates are in blue (solid). rates also matc h the true ones w ell. How ev er, after the a brupt c hange at time 12800s, the estimated flo w rates app ear to deviate from the true o nes. This p ossibly reflects the relativ e cum b ersomeness, to some exten t, of using constan t v ariances in the situations with rapid c hanges in the fluid dynamics. The deviation of the estimated flow rates from the true ones, in t erms o f the ro ot mean squared errors (RMSEs) av eraged o ve r the time in terv al b et w een 10000 s and 16000s, is 0.2975. 3.1.2 Results with the fixed interv al smo othing approac h T o minimize the appro ximate cost function in Eq. (6), w e adopt the MultiStart algorithm in the Global Optimization T o olb ox o f MA TLAB c  R2012a, with 3 start p o ints. The opti- mized v ar iance f or t he gas flow rate is 0.004 , while that for the o il flow rate is 0.1378. Fig. 5 sho ws the estimated flow ra tes (pa nel (a)) and the corresponding sim ulat ed temp erature and pressure (panel (b)), when the APF is used in conjunction with the flo w rate v ariances optimized through the fixed in terv al smo othing a pproac h. With the constan t estimated v ariances, the assimilation results in Fig. 5 a pp ear similar t o those in Fig . 4. How ev er, its corresp onding time mean RMSE becomes 0.2216, low er than t hat in Fig. 4. 3.1.3 Results with the lag-1 smo ot hing approac h No w w e examine the p erforma nce o f the APF when it is used in conjunction with the lag-1 smo othing approac h. T o implemen t the iteration pro cess of the EM alg o rithm, the initial 13 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 0.5 1 1.5 2 2.5 3 3.5 Estimated parameter p t , p a ± 2 σ Time [s] Zonal flow rate [kg/s] 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 9 10 11 12 13 14 15 16 17 Estimated parameter p t , p a ± 2 σ Time [s] Zonal flow rate [kg/s] (a) Flow rates 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3475 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 320 330 340 350 360 Sensor @ 3475 m, d o , d a Temperature [K] Time [s] 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3925 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 338 340 342 344 346 348 Sensor @ 3925 m, d o , d a Temperature [K] Time [s] (b) Pr e s sure and temp era ture Figure 5: As in Fig. 4, but the APF is used in conjunction with the flow rate v ar ia nces estimated by t he fixed in terv al smo othing approa c h. No t e that in this case the spread of the ga s flo w rate estimates is v ery small a nd is nearly in visible in the left figure of panel (a). v ariances of b oth the gas and oil flo w rates are set to 1. The stopping condition is either (a) the relativ e Euclide an norm c hange is less than 10 − 3 , or (b) the iteration num ber reaches 100. In the exp erimen ts, it is found that the iteration pro cess typically stops after 3–4 iterations, meaning that the v ariances estimate conv erges quite fast. F or illustration, in Fig. 6 w e rep ort the (final) estimated v ariances of the gas (left) and oil (righ t) flo w rates as functions of time. As o ne can see there, the estimated v ariances oscillate with t ime frequen tly , and ma y cross different orders of magnitudes in a short time. O ne consequence of suc h oscillations is that the cor r espo nding estimated flo w rates ma y app ear less smo oth in comparison with the fixed in terv al smo othing approac h, a s will be sho wn later. F or the APF equipp ed with the lag- 1 smo othing appro a c h, its estimated flo w rates a nd the corresp onding sim ulated temp erature a nd pressure are rep orted in Figs. 7 ( a) and 7(b), resp ectiv ely . Comparing Fig s. 5(b) and 7(b), one can see the APF with b oth smo othing approac hes yield close p erfor mance in terms of matching the temp erat ur e and pressure of the g a uges. On the other hand, the APF with t he lag-1 smo othing approach exhibits differen t b eha vior fr om that with the fixed interv al smo othing approach, in terms of the estimated flow rates. Because the v ariances in the lag- 1 smo othing appro ac h oscillate with time, the estimated flo w rates tend to appear less smo oth than those in the fixed in terv al smo othing approach. F or the same reason, though, af ter the abrupt c hang e at time 12800 s, the APF with the lag- 1 smo o thing approac h exhibits b etter abilit y to track the true flo w rates, a s can b e seen in Fig. 7(a). Ov erall, wh en the APF is equipp ed with the la g-1 smo othing approach, the corresp onding time mean RMSE b ecomes 0.1916, further low er than tha t in the fixed in terv al approach. 14 1 1.2 1.4 1.6 x 10 4 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 10 1 Time[s] Estimated variances of gas flow rates 1 1.2 1.4 1.6 x 10 4 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 10 1 10 2 10 3 Time[s] Estimated variances of oil flow rates Figure 6: Time series o f t he estimated v ariances of the gas (left) and oil (righ t) flow rates. 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 0.5 1 1.5 2 2.5 3 3.5 Estimated parameter p t , p a ± 2 σ Time [s] Q 1 [kg/s] 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 9 10 11 12 13 14 15 16 17 Estimated parameter p t , p a ± 2 σ Time [s] Q 2 [kg/s] (a) Flow rates 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3475 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 320 330 340 350 360 Sensor @ 3475 m, d o , d a Temperature [K] Time [s] 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3925 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 338 340 342 344 346 348 Sensor @ 3925 m, d o , d a Temperature [K] Time [s] (b) Pr e s sure and temp era ture Figure 7: As in F ig. 4, but no w the APF is used in conjunction with the flo w rate v ariances estimated by the lag-1 smoo t hing approac h. 15 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3475 m Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 320 330 340 350 360 Sensor @ 3475 m Temperature [K] Time [s] 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3925 m Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 338 340 342 344 346 348 Sensor @ 3925 m Temperature [K] Time [s] Figure 8: Measured pressure and temp erature from the ga uges, generated b y the well flow mo del with the finer resolution (5m). 3.2 The imp erfect observ ation scenario (IOS) In a recen t w or k (Lo r entzen et al., 2014b), the framew ork prop osed in Loren tzen et al. (2014a) is a pplied to study tw o f ull- scale field cases, in whic h v arious uncertain ties (e.g., observ at io n errors) in data assimilation are inv olv ed. F or our purp ose here, we consider a simplified IOS in whic h certain imp erfection may b e in tro duced to the syn thetic observ ation system (i.e., the well flow mo del) in order to ev aluate the p erformance of the APF with the v ariances obtained through differen t w a ys. The exp erimental settings in this section are the same as those in the POS, except for the p oten tial imp erfection in the observ ation system resulting f r o m the follo wing tw o sources of unce rta inties. One is tha t the w ell flow mo del used to generate the syn thetic observ a t io ns has a r esolution different from that of the we ll flow mo del adopted in da ta assimilation. Sp ecifically , to generate the syn thetic observ at io ns, the spatial resolution of the w ell flow mo del is 5m ( in md) along the w ell tra jectory . While in data assimilation, the spatial r esolution of the w ell flo w mo del is 50m instead. F o r comparison, in Fig. 8 w e sho w the measuremen ts generated b y the w ell flo w mo del with the finer resolution (the true gas a nd oil flow rates are the same as t hose in Fig. 3(a)), whic h are used in all our exp erimen ts b elow. C omparing Figs. 3(b) and 8, there are some clear differences in the t emp era t ure data. The other source o f uncertain ties is that the p df of the observ ation error v k in Eq. (1 ) ma y b e mis-sp ecified. Concretely , let W k b e the diag onal cov ariance matrix that is used to generate the syn thetic observ ations with the finer resolution (5m) in the w ell flow mo del (here W k is constructed in the same w a y a s in the POS). In data assimilation w e consider three cases, with the observ ation error cov ariance matrix b eing 0 . 5 × W k , 1 × W k and 2 × W k , respectiv ely . The refore in the cases with 0 . 5 × W k and 2 × W k , the observ a t ion error co v ariance matrices are mis-sp ecified. 16 T a ble 1: Time mean RMSE of the APF when the observ ation error co v ariance matrix in assimilation is 0 . 5 × W k , 1 × W k and 2 × W k , resp ectiv ely . man ual fixed interv al lag-1 Time mean RMSE (0 . 5 × W k ) 1.5026 1.4869 1.1885 Time mean RMSE (1 × W k ) 1.5228 0.9321 1.3533 Time mean RMSE (2 × W k ) 1.8201 1.2431 1.5020 Fig. 9 sho ws the p erformance of the APF when the observ ation error co v ariance matrix is 1 × W k . The upp er panels rep ort the estimated flow rates ( left) and the corresp o nding sim ulated observ ations (right), resp ectiv ely , when the APF uses the man ually chose n v ari- ances (0.5 for b oth gas and flow rates) in Loren tzen et al. ( 2 014a). Similarly , the middle and lo w er pa nels presen t the results of t he APF when the fixed-lag and lag-1 smo othing ap- proac hes a re adopt ed to optimize the v ariances. As one can see in F ig . 9, the APF can still matc h the pressu re data w ell in all cases. Ho w ev er, it s matching of the temperature data is w or sened compare to the cases in the POS. Accordingly , the corresp onding estimated flo w rates also exhibit relative ly large estimation errors. Indeed, with the manu ally c hosen v ariances, the time mean RMSE b ecomes 1.5228, while those with the v ariances opti- mized by the fixed-lag and lag- 1 smo othing approache s ar e 0.9321 and 1.35 33, resp ectiv ely . Therefore, in t erms of estimation accuracies, the APF with man ually c hosen v ariances under-p erforms those with the v ariances optimized b y b oth smo othing approac hes, and in this pa rticular case the fixed-lag smo othing approach leads to the low es t time mean R MSE. F or brevit y , the experiment results in terms of time mean RMSEs, with the observ ation error co v ariance ma t r ix b eing 0 . 5 × W k and 2 × W k , resp ectiv ely , are summarized in T able 1. As one can see, in b oth cases, the APF with ma nually chose n v ariances again under- p erforms the one with the v ariances optimized b y the tw o smo othing approach es. In the case with 0 . 5 × W k , the lag-1 smo othing approac h seems to result in the b est estimation, while in the case with 2 × W k , it is the fixed-lag smo othing approac h that a pp ears to ha v e the low est time mean RMSE. Compared to the results in the POS (cf. Figs. 4 – 7), the p erformance of the APF deteriorates in the IOS fo r all three observ ation error co v ariance ma t rices used in data assimilation. The reason of the p erformance downgrade may b e largely attributed to the imp erfection in the observ a tion system, since in the Bay esian estimation framew ork, it is essen tial for one to kno w the statistical distributions of bo th mo del and observ ation errors (see, for example, Eqs. (A.2) and ( A.3 )), while the imp erfection in t he observ a t io n system ma y mak e the estimation b ecome sub-optimal. In suc h circumstances, it may b e desirable for one to adopt robust da t a assimilation algorithms (see, for example, Luo and Hoteit, 2014, 2011; Simon, 2006) that are less sensitiv e to t he precise kno wledge of the mo del and/or observ ation error distributions. An inv estigation in this asp ect is, how ev er, b ey ond the scop e of the curren t w ork. 17 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 1 1.5 2 2.5 3 3.5 Estimated parameter p t , p a ± 2 σ Time [s] Q 1 [kg/s] 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 2 4 6 8 10 12 14 16 18 20 Estimated parameter p t , p a ± 2 σ Time [s] Q 2 [kg/s] (a) Flow rates (manual) 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3475 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 320 330 340 350 360 Sensor @ 3475 m, d o , d a Temperature [K] Time [s] 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3925 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 338 340 342 344 346 348 Sensor @ 3925 m, d o , d a Temperature [K] Time [s] (b) P ressure and temp era ture (man- ual) 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 0.5 1 1.5 2 2.5 3 3.5 Estimated parameter p t , p a ± 2 σ Time [s] Q 1 [kg/s] 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 0 5 10 15 20 25 Estimated parameter p t , p a ± 2 σ Time [s] Q 2 [kg/s] (c) Flow rates (lag -1) 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3475 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 320 330 340 350 360 Sensor @ 3475 m, d o , d a Temperature [K] Time [s] 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3925 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 338 340 342 344 346 348 Sensor @ 3925 m, d o , d a Temperature [K] Time [s] (d) Pressure and temp erature (la g-1) 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 0.5 1 1.5 2 2.5 3 3.5 4 Estimated parameter p t , p a ± 2 σ Time [s] Q 1 [kg/s] 1.1 1.2 1.3 1.4 1.5 1.6 x 10 4 0 2 4 6 8 10 12 14 16 18 Estimated parameter p t , p a ± 2 σ Time [s] Q 2 [kg/s] (e) Flow rates (fixed-la g) 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3475 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 320 330 340 350 360 Sensor @ 3475 m, d o , d a Temperature [K] Time [s] 1 1.2 1.4 1.6 x 10 4 1.2 1.4 1.6 1.8 2 x 10 7 Sensor @ 3925 m, d o , d a Pressure [Pa] Time [s] 1 1.2 1.4 1.6 x 10 4 338 340 342 344 346 348 Sensor @ 3925 m, d o , d a Temperature [K] Time [s] (f ) Pressure and temp erature (fixed- lag) Figure 9: Upp er panels (a,b): P erformance of the APF with the man ually c hosen v a riances (0.5 for b oth the oil and ga s flow rates). Middle panels (c,d): Performance o f the APF with the v ariances ( 0 .1270 for the gas flow rate and 0.1075 for the oil flo w rate) optimized b y the fixed-lag smo othing approac h. Low er panels (e,f ): P erformance o f the APF with the v ariances optimized b y the lag-1 smo othing approac h. 18 4 Conclus ion In this w ork w e studied the soft-sensing problem ba sed on a previously es tablished frame- w ork. The salient features of this framew or k include the transien t, instead o f steady- state, w ell flo w mo del, a Mark ov jump pro cess that has b etter abilit y in capturing abrupt c hanges in the fluid dynamics, and a mo dern estimation a pproac h, the a uxiliary particle filter (APF), that guides the estimated flo w rates tow ard matc hing the temp erature and pressure from the ph ysical sensors . This w ork fo cused on estimating the v ariances of the flo w rates in the Marko v jump pro cess. T o this end, w e prop o sed and implemen ted t w o approache s. One utilizes all the a v ailable measuremen ts, e.g., temp erature and pressure, in a fixed t ime in terv al, and is th us called the fixed interv al smo ot hing a ppro ac h. The other uses the measuremen ts up to (and including) 1-step a head of the time instan t at whic h the v ariances of the flow rates are to b e estimated, and therefore is named the lag-1 smo othing approach. In our implemen tations, the fixed in terv al smo othing approac h prov ides constan t v ariance estimates to the APF, and ma y th us tend to pro duce mor e smo oth flow rate estimation than the lag-1 smo othing approac h. On the other hand, the APF with the lag-1 smo othing appro ac h requires less w aiting time for the neces sary measure ments to arriv e, and tends to run faster in terms of computational time. In a nume rical example, w e inv estigated t he p erformance of our fra mew ork in b oth the p erfect o bserv ation scenario (POS) and the imp erfect o bserv ation scenario ( IOS) , in whic h the APF ra n with the man ually chosen v ariances, and those supplied b y the fixed in terv al and lag-1 smo othing a pproac hes. In t he POS, we sho we d that when t he APF is used in conjunction with the v ariances fro m either of the three metho ds, the sim ulated temp erature and pressure of the estimated flo w rates exhibit go o d matc hing to those of the ph ysical sensors. In terms of estimation accuracies o f the flo w rates, how ev er, the APF with the v ariances optimized by either smo o thing approac h seems to outp erform the one with manually c hosen v ariances. Similar results w ere also observ ed in the IO S. Ho w ev er, due to the imperfection in the observ ation sys tem, the p erformance of the APF tends to deteriora te, in t erms o f b oth t he estimation accuracy and the ability to match t he observ at io ns. T hese results thu s suggest the imp ort a nce f o r one to ha ve an observ atio n system tha t is able to model the underlying ph ysics reasonably w ell. Ac kno wled gmen t The authors w ould lik e to t ha nk tw o anony mous review ers for their constructiv e commen ts and suggestions that ha v e significan tly impro v ed the work. W e also ac know ledge financial supp orts fro m the Researc h Council o f Norwa y , Cono coPhillips Sk andina via AS a nd GDF SUEZ E&P Norge AS, through t he pro ject “T ransien t w ell flo w mo delling and mo dern estimation tec hniques for accurate pro duction allo cat io n.” 19 A The auxiliary particle filter Here w e provide an outline of the main pro cedures in the auxiliary particle filter (APF), follo wing Luo and Hoteit (2 014); Pitt and Shephard (1999). F or illustration, consid er the state estimation problem in the follow ing sy stem Q k = M k ,k − 1 ( Q k − 1 ) + u k , (A.1a) y k = H k ( Q k ) + v k . (A.1b) Here, Q k ∈ R r is the r -dimensional mo del state at time instan t k , y k ∈ R p the corre- sp onding observ a tion of Q k , u k ∈ R r the mo del error with the asso ciated p df p ( u k ), and v k ∈ R p the observ ation error with the asso ciated p df p ( v k ). The transition op erator M k ,k − 1 : R r → R r maps Q k − 1 to Q k (and it can b e a sto c hastic map as in Eq. (2)), and the observ ation op erator H k : R r → R p pro jects Q k from the state space on to the observ a- tion space. The problem of our in terest is to estimate the p osterior p df of the mo del state Q k at time instan t k , giv en the o bserv ations Y k = { y k , y k − 1 , · · · } up to and including k , together with the prior p df p ( Q i | Y i − 1 ) of the mo del state Q i at some earlier instan t i ( i ≤ k ). Recursiv e Ba y esian estimation (RBE) ( Arulampalam et al., 2002) provide s a proba- bilistic framew ork that recursiv ely solv es the state estimation problem in terms of some conditional p df s. Let p ( Q k | Y k − 1 ) b e the prior pdf of Q k conditioned on the observ atio ns Y k − 1 up to and including time k − 1, but without the know ledge of the observ a tion y k y et. Once the observ ation y k is kno wn, one incorpo rates the infor ma t ion conten t of y k ac- cording to Ba y es’ rule to up da t e the prior p df t o the po sterior one p ( Q k | Y k ). By ev olving p ( Q k | Y k ) forw ard in time, one obtains a prior p df p ( Q k +1 | Y k ) at the next time instan t. Concretely , the mathematical description of RBE consists of (Arulampalam et al., 2002): Prediction step: p ( Q k | Y k − 1 ) = Z p ( Q k | Q k − 1 ) p ( Q k − 1 | Y k − 1 ) d Q k − 1 , (A.2) and filtering step: p ( Q k | Y k ) = p ( y k | Q k ) p ( Q k | Y k − 1 ) R p ( y k | Q k ) p ( Q k | Y k − 1 ) d Q k , (A.3) where the transition p df p ( Q k | Q k − 1 ) a nd the like liho o d function p ( y k | Q k ) are assumed kno wn, in ligh t of the kno wledge of the pdfs o f the mo del and observ ation errors in Eq. (A.1). Once the explicit forms o f the conditional p dfs in Eqs. (A.2) a nd (A.3) are obtained, the o ptimal estimate and other asso ciated statistical information can b e derive d based on a certain optimalit y criterion, e.g., minimu m v ariance or maxim um lik eliho o d. Th us RBE pro vides a solution of the estimation problem, and conceptually leads to the optimal no nlinear filter. In practice, ho w ev er, difficulties often arise in deriving the exact optimal filter, largely due t o the fact that the in tegrals in Eqs. (A.2) and (A.3) are of ten intractable. Therefore 20 one may hav e to adopt a certain approxim at ion sc heme for ev aluation. In many situations, it if often reasonable to approx imate t he prior and p osterior p dfs b y certain Gaussian distributions. With this approxim at io n, the RBE fra mew ork reduces to the KF or certain extensions in nonlinear systems , e.g., the extended Kalman filter (EKF , see Jazwinski, 1970; Simon , 200 6 ), the sigma p oin t Kalman filters (SPKF, see, for example, Ito and Xiong , 2000; Julier et al., 2000; Luo and Moroz, 2009; Luo et al., 2012; Nørgaard et al., 2000), and the ensem ble Kalman filter (EnKF, see, for ex ample, Anderson, 2 001; Bishop et al., 2001; Burgers et al., 1 998; Luo and Hoteit, 2012, 2013; Pham, 2 001; Whitak er and Hamill, 2002). In addition, a family of nonlinear and no n- Gaussian data assimilation algorithms ba sed on the Gaussian mixture mo del (GMM) approx imatio n, often called the Ga ussian sum filter ( see, for example, Alspac h and Sorenson, 1972; Chen and Liu , 20 00; Hoteit et al., 2012, 20 08; Luo et al., 2010; Sorenson and Alspac h , 1971; Stordal et al., 2010), can b e constructed as a bank of parallel nonlinear Kalman filters (e.g., the EKF, the SPKF, or the EnFK), whic h enhances the re-usabilit y of existing codes of nonlinear Kalman filters. Sp ecific to our fo cus in the curren t w ork, if the Mon te Carlo approximation is adopted to approxim ate the prior and p o sterior p dfs in the R BE framew ork, then it leads to the par - ticle filter (PF, see, for example, Doucet et al., 2001; Gordon et al., 1993; Luo and Hoteit, 2014; Pitt and Shephard, 19 99) that is also applicable to no nlinear and non-G aussian data assimilation problems. F or instance, the p osterior p ( Q k − 1 | Y k − 1 ) a t the ( k − 1) t h step can b e approximated by p ( Q k − 1 | Y k − 1 ) ≈ n X i =1 w k − 1 ,i δ ( Q k − 1 − Q a k − 1 ,i ) , where Q a k − 1 ,i ( i = 1 , 2 , · · · , n ) are the particles at the filtering step b efore a re-sampling algorithm (if necessary) is applied, w a k − 1 ,i are the asso ciated w eigh ts, and n is the tota l n um b er of particles. F or notatio nal conv enience , let ˜ Q a k − 1 ,i b e the particles g enerated by a re-sampling algorithm, and ˜ w a k − 1 ,i the corresponding w eights. Supp o se that w e start with some part icles Q b k ,i from the prior p df at time instan t k , together with their asso ciated w eigh ts w b k ,i . Then, in consistency with Eqs. (A.2) and (A.3), the APF has the follo wing steps. • Filtering step: With an incoming observ ation y k , the particles remain unchanged, i.e., Q a k ,i = Q b k ,i , while the asso ciated weigh ts a r e up da t ed according to Ba ye s’ r ule so that w a k ,i = w b k ,i p ( y k | Q b k ,i ) n P i =1 w b k ,i p ( y k | Q b k ,i ) , (A.4) where p ( y k | Q b k ,i ) is the probability that y k happ ens to b e the observ a tion with resp ect to Q b k ,i . • Re-sampling step: In high- dimensional mo dels, w eigh t collapse ma y happ en and this will, in g eneral, deteriorate the p erformance of the PF. T o mitiga t e t his problem, 21 it is customary to include a re-sampling step that singles out a sp ecified n um b er of particles from the p o ol { Q a k ,i } n i =1 based o n a certain rule. In the con v entional PF, the selection is conducted based on the w eights { w a k ,i } n i =1 asso ciated with { Q a k ,i } n i =1 . Ho w ev er, as po in ted o ut b y Pitt and Shephard (1999), particles selec ted in this w a y are not guaran teed to match the future observ ations. As a remedy , one ma y tak e in to accoun t the particles’ abilit y to mat c h future observ ations. T o this end, one can, fo r instance, compute the inte rmediate w eights µ k ,i ∝ w a k ,i p ( y k +1 |M k +1 ,k ( Q a k ,i )), where p ( y k +1 |M k +1 ,k ( Q a k ,i )) represen ts the lik eliho o d function v alue of the for ecast M k +1 ,k  Q a k ,i  of Q a k ,i at time instan t k + 1. Based on µ k ,i , one applies sequen tial imp ortance re-sampling (SIR, see, for example, Arulampalam et al., 2002) to the indices i = 1 , 2 , ...n and obtains a new set of indices i s , s = 1 , 2 , · · · , n , to pic k up samples from the set { Q a k ,i } n j =1 . The indexed samples { ˜ Q a k ,i s } n s =1 then has equal w eigh ts, i.e., their asso ciated w eights ˜ w a k ,i s = 1 /n for all i s , s = 1 , 2 , · · · , n . • Predic tio n step: The re-sampled particles { ˜ Q a k ,i s } n s =1 are propagated f orw ard in time through the dynamical mo del Eq. (A.1a ) , so as to o btain a set of part icles { Q b k +1 ,s } n s =1 at time instant k + 1. Since at the previous re-sampling step, the lik eliho o d function v alues with respect to the observ ation at time k + 1 are inv olv ed in the re-sampling step, o ne needs to adjust the w eigh ts of { Q b k +1 ,s } n s =1 to mak e them consisten t with the R BE framew ork. T o this end, the we ights of Q b k +1 ,s are computed according to w b k +1 ,s ∝ p ( y k +1 | Q b k +1 ,s ) /p ( y k +1 |M k +1 ,k ( ˜ Q a k ,i s )). With these computed, one can pro ceed to the filtering step as describ ed previously , and so o n. References Aanonsen, S., G. Nævdal, D. Oliv er, A. Reynolds, and B. V all` es, The ensem ble Kalman filter in reserv oir engineering: a review, SPE Journal , 14 (3), 393–412, 200 9. Alspac h, D. L. and H. W. Sorenson, Nonlinear Ba y esian estimation using G aussian sum appro ximation, IEEE T r ansactions o n Automatic Contr ol , 17 , 439 – 448, 1972. Anderson, J. L., An ensem ble adjustmen t Kalman filt er for data assimilation, Mon. We a. R ev. , 129 , 288 4–2903, 2001. Arulampalam, M., S. Mask ell, N. G ordon, and T. Clapp, A tuto rial on particle filters f or online nonlinear/non- G aussian Ba y esian trac king, IEEE T r ans a ctions on Signal Pr o c ess- ing , 50 , 17 4 –188, 2002. Azim, R. R. A., Nov el applications of distributed t emp era t ur e measuremen ts to estimate zonal flow rate and pressure in offshore gas w ells, in The SPE International Pr o duction and Op er ations Confer enc e an d Exhibition , Doha, Qatar, 2 012, SPE 154098. 22 Bishop, C. H., B. J. Etherton, and S. J. Ma jumdar, Adaptiv e sampling with ensem ble transform Kalman filter. Part I: theoretical asp ects, Mon. We a. R ev. , 129 , 420 – 436, 2001. Blo emen, H. H. J., S. P . C. Belfroid, W. L. Sturm, and F. J. P . C. M. G. V erhelst, Soft sensing for gas-lift w ells, SPE Journal , 11 , 454 – 463, 2006, SPE 90370. Burgers, G., P . J. v an Leeu w en, and G. Ev ensen, On the analysis sc heme in the ens emb le Kalman filter, Mon. We a. R ev. , 126 , 1719–1724, 1998. Chen, R. a nd J. Liu, Mixture Ka lma n filters, Journal of the R o yal Statistic al S o ciety: Series B (Statistic al Metho dolo gy) , 62 (3), 49 3–508, 2000. de Kruif, B., M. Lesk ens, R. v an der Linden, and G. Alberts, Soft-sensing for multilateral w ells with do wn hole pressure and temp erature and surface flo w measureme nts, in A bu Dhabi I nternational Petr oleum Exhibition and Confer enc e , Abu Dhabi, UAE, 20 0 8, SPE 118171. Dempster, A., N. Laird, and D. Rubin, Maxim um lik eliho o d from incomplete data via the EM algorithm, Journal of the R oyal Statistic al So ciety, S e ri e s B , 39 , 1–3 8, 1977. Doucet, A., N. D e F reitas, and N. Gordon (eds.), Se quential Monte Ca rlo metho ds in pr actic e , Springer V erlag, 2001. Ev ensen, G., Data Assimilation: The Ensemble Kalman Filter , Springer, 2006. Gordon, N. J., D. J. Sa lmond, and A. F. M. Smith, No v el approa ch to nonlinear a nd non- Gaussian Bay esian state estimation, IEE Pr o c e e dings F in R adar and Sig nal Pr o c essing , 140 , 107–11 3, 1993. Gryzlo v, A., W. Schiferli, and R. F. Mudde, Estimation of multiphase flow rates in a horizon tal w ellb ore using an ensem ble K a lman filter, in 7th International Confer enc e on Multiphase Flow ICMF , T ampa, Flor ida , USA, 20 1 0. Hoteit, I., X. Luo, and D. T. Pham, P article Kalman filtering: A n optimal nonlinear framew ork f or ensem ble Kalman filters, Mon. We a. R ev. , 140 , 528–542, 201 2 . Hoteit, I., D. T. Pham, G. T r ia n tafyllou, and G. Korres, A new approxim at e solution o f the optimal nonlinear filter for dat a assimilation in meteorology and o ceanograph y , Mon. We a. R ev. , 136 , 317–334, 2008. Ito, K. and K. Xio ng , Gaussian filt ers f o r nonlinear filtering problems, IEEE T r ansactions on A utomatic Contr ol , 45 , 9 10–927, 200 0 . Jazwinski, A. H., S to chastic Pr o c esses and Filtering The ory , Academic Press, 1 970, 400 pp. 23 Julier, S., J. Uhlmann, and H. D urran t-Whyte, A new metho d f or the nonlinear t r a ns- formation of means and cov ariances in filters and estimators, IEEE T r ansactions on A utomatic Contr ol , 45 , 47 7–482, 2000 . Lesk ens, M., B. de Kruif, a nd S. B. J. Smeulers, D o wnhole m ultiphase metering in wells b y means o f soft-sensing, in Intel lige nt Ener gy Confer enc e and Exhibition , Amsterdam, Netherlands, 20 0 8, SPE 11 2 046. Loren tzen, R ., A. Stordal, G. Nævdal, H. Karlsen, and H. Sk aug, Estimation of pro duction rates using t ransien t w ell flow mo deling a nd the auxiliary part icle filt er, SPE Journal , 19 , 172 – 180, 2014a. Loren tzen, R. J., O. Sæv areid, and G. Nævdal, Rat e allo cation: Com bining transien t w ell flo w mo deling and data assimilation, in T he SPE A nnual and T e chnic al Confer enc e and Exhibition , Florence, Italy , 2010a, SPE 13 5073. Loren tzen, R. J., O. Sæv areid, and G. Nævdal, Soft m ultiphase flo w metering for accurate pro duction allo catio n, in T he SPE Russian O i l & Gas T e chnic al Confe r enc e , Mosco w, Russia, 2010 b, SPE 1360 2 6. Loren tzen, R. J., A. S. Stordal, X. Luo, and G. Nævdal, Estimation of pro duction rates using transient w ell flo w mo deling and the auxiliary pa rticle filter - full- scale a pplicatio ns, SPE Pr o duction & Op er ations , in review , 2014b. Luo, X. and I. Hoteit, Robust ensem ble filtering and its relation to co v ariance inflation in the ensem ble Kalman filter, Mon. We a. R ev. , 139 , 3938–3953, 2011. Luo, X. and I. Hoteit, Ensem ble Kalman filtering with residual nudging, T el lus A , 64 , 17,130, op en access , 2012. Luo, X. and I. Hoteit, Cov ariance inflation in the ensem ble Kalman filter: a residual n udging persp ectiv e and some implications, Mon. We a. R ev. , 141 , 33 60–3368, 20 1 3. Luo, X. and I. Hoteit, Efficien t particle filtering through residual n udging, Quart. J. R oy. Mete or. So c. , 140 , 557–572, 2014 . Luo, X., I. Hoteit, and I. M. Moroz, On a nonlinear Kalman filter with simplified divided difference approximation, Physic a D , 241 , 6 71–680, 201 2 . Luo, X. and I. M. Moro z, Ensem ble Kalman filter with the unscen ted transform, Physic a D , 238 , 54 9 –562, 20 09. Luo, X., I. M. Mor o z, and I. Hoteit, Scaled uns cen ted transform Gaussian su m filte r: Theory and application, Physic a D , 239 , 684 – 701, 201 0. Nævdal, G., L. M. Johnsen, S. I. Aanonsen, E. H. V efring, et al., Reserv oir monitoring a nd con tin uous mo del up dating using ensem ble K alman filter, SPE journal , 10 , 66–74 , 20 05. 24 Nørgaard, M., N. K. P oulsen, and O. Ravn, New dev elopmen ts in state estimation for nonlinear systems, A utomatic a , 36 , 1 627 – 1 638, 200 0. Pham, D. T., Sto chastic metho ds for sequen tial dat a assimilation in strongly nonlinear systems , Mon. We a. R ev. , 129 , 1194 –1207, 2001. Pitt, M. and N. Shephard, Filtering via sim ulation based a uxiliary particle filters, J. Amer- ic an S tatistic al Asso ciation , 94 , 590–599, 1999. Silv erman, B. W., Density Estimation for Statistics and Data Analysis , Chapman & Hall, 1986. Simon, D., Optimal State Estimation: Ka l m an, H-Infinity, and Nonline ar Appr o aches , Wiley-In terscience, 2006, 552 pp. Sorenson, H. W. and D. L. Alspac h, Recursiv e Ba y esian estimation using G a ussian sums, A utomatic a , 7 , 465 – 479, 1971. Stordal, A. S., H. A. Karlsen, G. Nævdal, H. J. Sk aug, a nd B. V all` es, Bridging t he ensem ble Kalman filter and pa r ticle filters: the adaptiv e Gaussian mixture filter, Computational Ge oscienc es , 15 , 293–305, 2010 . Whitak er, J. S. and T. M. Hamill, Ensem ble data assimilation without p erturb ed observ a- tions, Mon. We a. R ev. , 130 , 1913– 1924, 2002. W r ob el, K. a nd W. Schiferli, Soft-sensing, non-in trusiv e multiphas e flo w meter, in 14th International C onfer enc e o n Multiphase Pr o duction T e chnolo gy , Cannes, F ra nce, 2009, SPE 2009- B1. Y oshiok a, K., D. Zh u, A. Hill, and L. W. Lak e, A new inv ersion method to in terpret flo w profiles fro m distributed temp erature and pressure measuremen ts in horizon tal wells , SPE Pr o duction & Op er ations , 24 , 510–5 2 1, 2009. 25