We focus on the performance of the xed stress splitting algorithm for the immiscible water-oil ow coupled with linear poromechanics. The two-phase ow equations are solved on general hexahedral ele- ments using the multipoint ux mixed nite element method, whereas, the geomechanics equations are dis- cretized using the conforming Galerkin method. The e ects of the coupling parameter on the performance of the xed stress algorithm is demonstrated for the Frio oil reservoir site. Keywords Immiscible two phase water-oil ow  poromechanicscoupled system of equations  staggered solution algorithm mixed nite element methodconforming Galerkin method 1 Introduction The understanding of the physics of subsurface phe- nomena is critical to the study of a multitude of applica- tions in the environment and the energy industry such as CO 2sequestration [1], hydrocarbon extraction [2], water ooding [3] and hydraulic fracturing [4] among many others. This physics is typically governed by the Saumik Dana Mork Family Department of Chemical Engineering and Ma- terial Science, University of Southern California, CA 90007 Mohamad Jammoul Oden Institute for Computational Engineering and Sciences, The University of Texas at Austin, TX 78712 Mary F Wheeler Oden Institute for Computational Engineering and Sciences, The University of Texas at Austin, TX 78712theory of multiphase ow through porous media [5, 6]. The in silico resolution of such physics requires the so- lution of a complicated set of equations governing the response of the various parameters under the appli- cation of thermo-mechanical-hydraulic loading [7{39]. Staggered solution schemes are those in which an op- erator splitting strategy is used to split the coupled problem into well-posed subproblems which are then solved sequentially instead of solving the coupled sys- tem of equations monolithically. The xed stress split is one such strategy that solves the ow problem while freezing the mean stress followed by the poromechanics problem in repeated iterations until a certain conver- gence criterion is met at each time step [40{63]. In this paper, we present the xed stress formulation for the immiscible two-phase (water-oil) ow coupled with linear poromechanics. The xed stress splitting leads to two updates for the Lagrangian porosity, a uid update and a mechanics update. The convergence criterion for the xed stress splitting aims to reduce the di erence between the two quantities. We study the error patterns in the porosity update re- sulting from adopting a xed tolerance in the conver- gence criterion of the splitting algorithm. This document is structured as follows: the model equa- tions are presented in this Section, the xed stress split algorithm is presented in Section 2, a numerical result for a eld case showing the performance metrics is pre- sented in Section 3 and the conclusions and outlook are presented in Section 4. 1.1 Two-Phase Flow Model Let R3be the ow domain with boundary @ = f D[f Nwheref Dis Dirichlet boundary and f Nis 2 Saumik Dana et al. Neumann boundary. The equations of ow model are @(  S ) @t+rz @(  S ) @t=q (mass conservation ) z =K (rp  g) (Darcy0s law )  = 0ec (p p 0)(slightly compressible fluid ) pc=popw (capillary pressure ) Sw+So= 1 (constraint ) z n= 0 onf N(0;T] p (x;0) =p 0(x); (x;0) = 0(x)8x2 where wfor the water (wetting) phase, ofor the oil (non-wetting) phase, p : (0;T]!Rare the phase pressures, z : (0;T]!R3are the phase uxes,S are the phase saturations,  are the phase densities, = (1 +) is Lagrangian porosity where is the Eulerian porosity, and is the volumetric strain, nis the unit outward normal on f N,q is the source or sink term, Kis the uniformly symmetric positive def- inite absolute permeability tensor,  =k  = are the phase mobilities where k are the phase relative permeabilities and  are the phase viscosities, c are the uid compressibilities and T > 0 is the time inter- val. The reader is refered to [63] for an understanding of mixed nite element formulation of ow model. 1.2 Linear Poromechanics Model Let R3be the poromechanics domain with bound- ary@ =p D[p Nwherep Dis Dirichlet boundary and p Nis Neumann boundary. The equations are r+f=0 (linear momentum balance ) =D pI (constitutive law ) f=  g+r(1 )g p=pwSw+poSo=popcSw(average phase pressure ) =wSw+oSo (average phase density ) (u) =1 2(ru+rTu) (small strain kinematics ) un1= 0 onp D[0;T];Tn2=tonp N[0;T] u(x;0) =08x2 where u: [0;T]!R3is the solid displacement, ris the rock density, is the Biot parameter, fis body force per unit volume, tis the traction boundary condition, is the strain tensor, Dis the fourth order elasticity tensor and Iis the second order identity tensor. The mean stress vis given as v=1 3tr() =1 3tr I+ 2G pI yes noSolveﬂow with mean stress ﬁxedNew coupling iteration Solve poromechanics converged?New time stepFig. 1: Flowchart for xed-stress split iterative scheme for two-phase ow coupled with linear poromechanics. =1 33+1 32G * tr()1 33 p = +1 32G  p=Kb p whereKbis the bulk modulus of the solid skeleton or the drained bulk modulus. The reader is refered to [63] for an understanding of continuous galerkin formulation of poromechanics model. 2 Fixed Stress Split Algorithm The owchart of the algorithm is provided in Figure 1. The basic idea of the xed - stress split strategy is to solve for the pressure and ux degrees of freedom at the current xed - stress iteration based on the value of the mean stress from the previous xed - stress iteration. These pressures then contribute to the force vector in the poromechanics system which is solved for displace- ments thereby updating the stress state. This updated stress state is then fed back to the ow system for the next xed - stress iteration. The stopping criterion for coupling iterations is designed based on the require- ment that the change in Lagrangian porosity during the poromechanics solve in every coupling iteration re- duces monotonically from one coupling iteration to the next. Now the change in Lagrangian porosity during the poromechanics solve in any coupling iteration in every time step is given by  por  flow where  poris the Lagrangian porosity at the end of the poromechanics solve,  flow is the Lagrangian porosity at the end of the ow solve. This quantity is normalized by the cur- Performance metrics of the xed stress split algorithm for multiphase poromechanics 3 rent value of Lagrangian porosity as  por  flow porand the stopping criterion is based on the requirement that the maximum norm of this quantity is less than a certain pre-speci ed tolerance as follows  por  flow por 1TOL 2.1 Lagrangian Porosity Update During Flow Solve The nonlinear nature of the equations of two phase ow neccesitate the need for nonlinear iterations in the man- ner of the Newton-Raphson method. The xed mean stress constraint is imposed in every Newton iteration. The evolution law for Lagrangian porosity is given by (k) = (k)p (2.1) where(k)() refers to the change in quantity ( ) over the (k+ 1)thNewton iteration and is the rock com- pressibility given by = Kb+ (1 +)  Kb(2.2) The reader is refered to [63] for the derivation of (2.1). 2.2 Lagrangian Porosity Update During the Poromechanical Solve The evolution law for Lagrangian porosity is given by = 0+ +(  0)(1 ) Kb(pp0) (2.3) The reader is refered to [63] for the derivation of (2.3). 3 Numerical result In this section, we discuss the numerical results for the Frio oil eld site. The simulation was performed us- ing the in-house reservoir simulator IPARS ( https:// csm.ices.utexas.edu/IparsWeb/index.html ) on the Bevo3 cluster at the Oden Institute for Computational Engineering and Sciences at the University of Texas at Austin. 3.1 Frio Field Case The Frio site is an oil eld located on the Gulf Coast of the United States near Dayton, Texas. The eld has a challenging geometry with several thin faults cross- ing through. The whole eld is curved in the vertical direction with a depth di erential of more than 1,000 ft between its parts [64]. Y XZXPERM 1600 1400 1200 1000 800 600 400 200(a) x-permeability (mD) Y XZZPERM 180 160 140 120 100 80 60 40 20 (b) z-permeability (mD) Fig. 2: Example 1: Permeability in the x-direction (left) and z-direction (right) for the Frio eld. 3.1.1 Con guration and Parameters We used hexahedral elements to capture the complex reservoir geometry. We have also incorporated the avail- able petrophysical and mechanical information to build a 3D model inside IPARS. Figures 2 show the eld ge- ometry and the distribution of the permeability in the x-direction and z-direction relatively. The computational domain was discretized with 33 9 3 elements. No- ow boundary conditions were imposed for the uid equations and both stress and displace- ment boundary conditions were used for the mechanics equations. We assumed a constant porosity = 0:20 ev- erywhere in the reservoir. The input parameters for the problem are summarized in Table 1. Three production wells and six injection wells were considered in this ex- ample. The locations of the wells are shown in Figures and . All nine wells have constant bottomhole pressure; the values of which are included in Table 1. 4 Saumik Dana et al. Quantity Value Time step 0:05days Number of grids 891(33 93) Absolute permeabilities (5:271010to3:1103)md Initial porosity 0:2 Water viscosity 1:0cp Oil viscosity 2:0cp Initial oil concentration 44:8lbm=ft3 Initial oil pressure 400:0psi Water compressibility 1E6psi1 Oil compressibility 1E4psi1 Rock density 165:44lbm=ft3 Initial water density 56:0lbm=ft3 Initial oil density 62:4lbm=ft3 Young's Modulus 1:2E6psi Poisson's ratio 0.35 Flow boundary conditions No ow on all boundaries X+" boundary xx=nx= 10000psi X-" - boundary u=0, zero displacement Y+" - boundary u=0, zero displacement Y-" - boundary yy=ny= 2000psi Z+" - boundary u=0, zero displacement Z-" - boundary zz=nz= 1000psi Injection well (1): Pressure speci ed, 4000.0 psi Injection well (2): Pressure speci ed, 3300.0 psi Injection well (3): Pressure speci ed, 4000.0 psi Injection well (4): Pressure speci ed, 4400.0 psi Injection well (5): Pressure speci ed, 3700.0 psi Injection well (6): Pressure speci ed, 4400.0 psi Production well (1): Pressure speci ed, 2000.0 psi Production well (2): Pressure speci ed, 2000.0 psi Production well (3): Pressure speci ed, 2000.0 psi Table 1: Parameters for the Frio eld problem 3.1.2 Results The initial reservoir state and end of simulation results for the Frio eld are displayed in Figures 3 and 4 respec- tively. The Figures (4a - 4d) show the reference phase pressure, water saturation, oil concentration, and ver- tical displacement after 60 days. The areas near the wells show signi cant change in pressure, water satu- ration, and oil concentration (Figures 4a and 4b). The change in the pressure induced surface subsidence in the uppermost corner of the reservoir as shown in Fig- ure (4d). We attempt to study the performance of the split- ting algorithm by looking at the error at the beginning and end of each xed stress iteration (Figures 5a and 5b). The error at the rst xed stress iteration of each time step decreases gradually. This can be correlated to the change in pressure in the reservoir. Since all the in- jection and production wells are pressure speci ed, the change in the pressure from one time step to another decreases gradually. This leads to a decrease in the er- ror at the start of the xed stress iteration as long as PWAT 5000 4600 4200 3800 3400 3000 2600 2200 1800 1400 1000(a) Reference phase pressure (psi) SWAT 0.4 0.375 0.35 0.325 0.3 0.275 0.25 (b) Water saturation COIL 75 70 65 60 55 50 45 (c) Oil concentration (%) DISPZ 2 1.5 1 0.5 0 ­0.5 ­1 (d) Vertical displacements (in) Fig. 3: Simulation results for the Frio eld at time t= 0:1 days. The locations of the injection (blue) and production (red) wells are shown. Performance metrics of the xed stress split algorithm for multiphase poromechanics 5 PWAT 5000 4600 4200 3800 3400 3000 2600 2200 1800 1400 1000 (a) Reference phase pressure (psi) SWAT 0.4 0.375 0.35 0.325 0.3 0.275 0.25 (b) Water saturation COIL 75 70 65 60 55 50 45 (c) Oil concentration (%) DISPZ 2 1.5 1 0.5 0 ­0.5 ­1 (d) Vertical displacements (in) Fig. 4: Simulation results for the Frio eld at time t= 60 days. The locations of the injection (blue) and production (red) wells are shown.the boundary conditions and source functions of the mechanics equations are independent of time. The error at the convergence of the xed stress split- ting for each time step is plotted in Figure 5b. The error admits aV-shape cycle and then decreases monotoni- cally afterwards. To understand this behavior, we plot the number of xed stress iterations for each time step. The increase in error occurred after the number of xed stress iterations decreased to a single iteration. The er- ror after the rst xed stress iteration, at the 7thtime step and onwards, was lower than the TOL = 104 for all values, and thus a single iteration was su- cient for the convergence of the algorithm. These ob- servations are consistent with the results of [54] where they reported that using the relative change in mean stress (and equivalently porosity) as a stopping crite- rion for the xed-stress iterations may easily lead to over-iteration (or under-iteration), unless the conver- gence threshold is carefully tuned. Another aspect to consider is the e ect of the cou- pling parameter, Biot's , on the performance of the splitting algorithm. By comparing the ve plots in Fig- ure 5c, increasing the value of leads to an increase in the error in the splitting algorithm and to decrease in the computational eciency. The decrease in compu- tational eciency can be inferred from the additional xed stress iterations needed to reach convergence for higher values. This behavior is expected because the increase in Biot's constant means a stronger coupling between the uid ow and mechanics. This is evidenced by the fact that the rock compressibility measuring the rate of change of Lagrangian porosity with respect to pore pressure increases with increase in Biot constant as shown in Equation (2.2). This also implies that the porous rock is more sensitive to changes in uid pres- sure inside the pores as the Biot constant increases. 4 Conclusions and Outlook The model equations and algorithmic framework for immiscible two phase ow coupled with linear porome- chanics have been presented in this work. We also stud- ied the convergence rates for the xed stress split method and the impact of the stopping criterion on the perfor- mance of the algorithm. The impact of Biot's parame- ters, namely , on the convergence and computational eciency of the xed stress algorithm was also consid- ered in this work. 6 Saumik Dana et al. 0 1 2 3 4 5 6 7 8 9 10 Timestep Number10-610-510-410-310-210-1100 (a) Error after rst xed stress iteration 0 1 2 3 4 5 6 7 8 9 10 Timestep Number10-610-510-410-310-210-1100 (b) Error after last xed stress iteration 0 1 2 3 4 5 6 Fixed Stress Iteration10-510-410-310-210-1100 (c) Fixed stress convergence for the rst time step 0 1 2 3 4 5 6 7 8 9 10 Timestep Number0123456Number of Fixed Stress Iterations (d) Number of xed stress iterations Fig. 5: Fixed stress metrics for the Frio eld simulationsCon ict of interest The authors declare that they have no con ict of inter- est.