Integrated field operation and optimization

Journal of Petroleum Science and Engineering 81 (2012) 161–170

Contents lists available at SciVerse ScienceDirect

Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol

Integrated field operation and optimization Silvya Dewi Rahmawati a,⁎, Curtis Hays Whitson b, Bjarne Foss a, Arif Kuntadi b a b

Department of Cybernetics Engineering, Norwegian University of Science and Technology, NO-7034 Trondheim, Norway Department of Petroleum Engineering and Applied Geophysics, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway

a r t i c l e

i n f o

Article history: Received 23 March 2011 Accepted 24 December 2011 Available online 2 January 2012 Keywords: petroleum field asset integrated optimization model non-derivative optimization method

a b s t r a c t Integrated modeling and optimization is an important method in the petroleum industry, particularly for field development and continuous asset-management evaluation. Traditional modeling consists of the application of many independent models describing the various elements of a petroleum asset in a silo-model approach, such as reservoir models, well models, surface process models, export and sales models and economic models. The effort to integrate these models such that the overall system performance can be optimized presents many technological challenges. A common field asset is modeled with emphasis on: (1) Three different reservoirs (two gas condensate reservoirs and one miscible oil reservoir), (2) A surface facility model (focusing on water handling, NGL extraction, sales-gas spec, and gas reinjection), and (3) Economic model. The surface process model interacts with the three reservoir models through distribution of available produced gas for reinjection into the three reservoirs. The multi-field asset model is introduced as a benchmark case. A decision support tool is used to develop the model and provides long-term production forecasts of gas, oil, and NGL revenue. Cost functions are introduced for all major control variables (number of wells, surface facility selection and operating conditions, and injection gas composition). Net present value is used as the target objective function. This paper evaluates optimal production strategies using several key control variables and field operational constraints. Optimization performance is tested with a non-derivative optimization solver, a constrained implementation of Nelder–Mead Simplex Reflection method. The benchmark is provided to the industry through a website containing application data files, network infrastructure, and results from our integrated optimization model. © 2011 Elsevier B.V. All rights reserved.

1. Introduction Operation of complex assets may require a holistic view of the value chain. This is particularly important if the different parts of the value chain are tightly connected. Current industrial practice typically takes a silo approach in the sense that one part of the supply chain is treated quite separate from other parts. This is particularly the situation in the upstream domain where, for instance, a decision support application for optimally allocating well production may include well and pipeline models. The downstream boundary condition is typically a constant pressure at the inlet separator. Similarly an optimizer for the surface process does not include models of the upstream system. This implies that the inlet separator acts as a “dividing wall” between two optimizers even though the two subsystems might be tightly connected. An example of this is when the gas output from the surface facility is fed back into the upstream system through gas-lift wells or gas injectors. There are many reasons for the silo-like approach currently ⁎ Corresponding author. Tel.: + 47 96881572. E-mail addresses: [email protected], [email protected] (S.D. Rahmawati). 0920-4105/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.petrol.2011.12.027

used in the industry. Different parts of the supply chain recruit people with different backgrounds and they use quite different decision support tools. This limits integration even in situations where integration has an obvious potential. Several researchers have conducted research on various integration topics. An integrated application for reservoir-production strategies and field development management has been developed by coupling an ECLIPSE® reservoir simulator with a surface and production network simulator and optimizer (Netopt) (Hepguler and Barua, 1997; Hepguler et al., 1997). (Trick, 1998) applied a somewhat different procedure from (Hepguler et al., 1997) using the same interface. In this case an ECLIPSE black oil reservoir simulator was coupled to a surface gas deliverability forecasting model, FORGAS. A discussion on production-injection operation systems (PIOS) based on concepts of system engineering has been proposed (Yang et al., 2003). (Lasschuit and Thijssen, 2004) discussed about supply chain planning and scheduling in the oil and chemical industry employing a mixedinteger non-linear programming (MINLP) related to costs, pricing and cargo costs. Cullick et al. (2003) and Bailey et al. (2005) discussed complex petroleum field projects applying uncertainty analysis, without considering a surface process model.

162

S.D. Rahmawati et al. / Journal of Petroleum Science and Engineering 81 (2012) 161–170

Tomasgard et al. (2007) presented a natural gas value chain model and integration applying an upstream perspective and a stochastic portfolio optimization. Run time can be a challenge in integrated applications, especially when closely linked high-fidelity models are tightly connected. The use of integrated optimization in day-to-day operations of the LNG value chain was studied by Foss and Halvorsen (2009). To reduce computation time they chose simple models for all system components. A sizable gain could be identified by integrating all models into one decision support application as opposed to dividing them into two applications; one for the upstream part and the other for the LNG plant. The literature cited above identifies a potential for integrating models in decision support tools. Moreover, integrated simulation and optimization is clearly regarded as an interesting but challenging topic. Hence, in this paper we present a benchmark problem which is designed to assess the potential of an integrated approach in decision support tools. A realistic benchmark as well as a base case will be defined in the following sections. Further, a sensitivity analysis of key decision variables will be presented in addition to some optimization results. 2. Integrated model The model presented in this paper is rich and complex enough to represent the value chain from reservoir to market and thus suitable as a benchmark for integrated operations and optimization (I-OPT). The upstream part of the I-OPT model includes two gas-condensate reservoirs and an oil reservoir while the surface process system includes gas and liquid separation as well as an NGL plant. The model also includes an economic component as indicated in Fig. 1. The project is designed with close links between the upstream and downstream parts of the model, partly due to gas re-injection. This is important since the I-OPT model is designed to study and assess the business value of integrated optimization as a decision support method. Integrated optimization in this context is defined as applications which utilize several different models along the value chain, for instance a reservoir model and a surface process model, in one optimization-based application as opposed to two separate applications for the reservoir and surface parts, respectively. Hence, the I-

OPT model is designed to challenge the conventional silo approach. The I-OPT model is further designed to study decisions both on a life-cycle horizon as well as shorter time frames. The surface facility model is steady-state while dynamic reservoir models account for depletion and displacement. The I-OPT model will be presented in the following sections. Complete documentation of the I-OPT model, including the base case discussed later, are made available at http:// www.ipt.ntnu.no/~io-opt/wiki/doku.php. 2.1. Reservoir model The reservoir models include two gas-condensate reservoirs and an oil reservoir. The gas-condensate reservoirs are scaled up from Juell et al. (2010) and the oil reservoir is a scaled up version of a miscible WAG project (Killough and Kossack, 1987). In the base case each reservoir is producing through 5 production wells and injection is made with 8 injection wells which provide gas cycling in the gascondensate reservoirs and miscible WAG in the oil reservoir. The production and injection wells are perforated through all layers. The size of the gas-condensate reservoirs is 3218 × 3218 × 48.7 m 3 which is divided into 36 × 36 × 4 grid blocks. The oil reservoir size is 5334 × 5334 × 30.48 m 3, divided into 35× 35× 3 grid blocks. The horizontal permeability distributions for the three reservoirs vary from a low permeability in the south west region towards higher permeability in the north east. The reservoir profiles are shown in Fig. 2 and the permeability distribution range is presented on Table 1. There are two faults in the horizontal direction, one is sealing and the other is partially communicating. The non-communicating fault separates low permeability and medium permeability areas. The partially communicating fault separates the medium and high permeability areas. The non-communicating shale in the vertical direction occurs between layers 3 and 4 in the lean gas-condensate reservoir, between layers 1 and 2 in the rich gas-condensate reservoir and between layers 2 and 3 in the oil reservoir. The reservoir models are compositional. The composition for the gas-condensate reservoirs consists of 9 components and the composition for the oil reservoir consists of 6 components. The initial fluid composition for each reservoir is presented in Table 2. The compositional reservoir models are run using the SENSOR® reservoir simulator.

IO Collaboration and Visualization Center Reservoir D e c i s i o n s u p p o r t

Reservoir & Well Simulator (Eclipse, VIP, Sensor)

Model

Pipelines

Network Simulator (GAP, OLGA, Wellflo)

Model

Surface Process

Surface Simulator (HYSYS, UniSim)

Model

Economic

NPV Calculation

Model

Integrated Operations and Optimization (I-OPT) Fig. 1. Integrated optimization schematic.

S.D. Rahmawati et al. / Journal of Petroleum Science and Engineering 81 (2012) 161–170

Lean Gas Condensate

a) Lean Gas Condensate Reservoir.

Table 1 Horizontal permeability and thickness distributions.

GINJ 2

GINJ 8

GINJ 4

Layer

PROD 3

PROD 2 GINJ 5

GINJ 1

PROD 5 GINJ 6

PROD 1

PROD 4

GINJ 3

GINJ 7

163

Sealing fault

1 2 3 4

Permeability (md)

Thickness (m)

Lean gas condensate

Rich gas condensate

Oil

Gas condensate

Oil

13–1300 4–400 2–200 15–1500

35–3500 4.5–450 2.5–250 1–100

50–5000 5–500 20–2000 –

9.1 9.1 15.2 15.2

6.1 9.1 15.2 –

TransX 0.00

20.46

40.92

61.37

81.83

Rich Gas Condensate

b) Rich Gas Condensate Reservoir.

Component

GINJ 5

PROD 2

GINJ 4

PROD 3

GINJ 1

GINJ 6 PROD 5 PROD 1 PROD 4

GINJ 7

Sealing fault

GINJ 3

28.64

57.28

85.92

114.56

c) Oil Reservoir. GINJ 2

GINJ 1 PROD 5

PROD 3

PROD 2

GINJ 8

Oil

Lean GC reservoir

GINJ 2 GINJ 8

TransX 0.00

Table 2 Initial composition for each reservoir.

PROD 1

GINJ 5

GINJ 6 PROD 4

GINJ 4

GINJ 3 Sealing fault

GINJ 7

CO2 0.01195 N2 0.01995 C1 0.66936 C2 0.10868 C3 0.06474 C4-6 0.07976 C7P1 0.03272 C7P2 0.01052 C7P3 0.00234 Oil reservoir (component) C1 C3 0.5 0.03

Rich GC reservoir Layer 1

Layer 2

Layer 3

Layer 4

0.05794 0.01852 0.62858 0.08288 0.0564 0.09251 0.04521 0.01468 0.00327

0.05789 0.01829 0.62398 0.0831 0.05691 0.09382 0.0469 0.01543 0.00367

0.05777 0.01789 0.61573 0.08347 0.05778 0.09611 0.04995 0.0168 0.00449

0.05713 0.01655 0.58684 0.08442 0.06048 0.10359 0.06078 0.02175 0.00845

C6 0.07

C10 0.2

C15 0.15

C20 0.05

During a simulation, the minimum tubing head pressure (THP) for each production well becomes a constraint. The THP for each well is compared with the manifold pressure from the surface calculation and is redefined as THP. The reason for this is to change the minimum THP to equal the manifold pressure when the manifold pressure is greater than the THP. The manifolds refer to the incoming streams in Fig. 3 (left side) for each reservoir. The distance between wellhead and manifold is assumed small therefore the pressure drop could be neglected. 2.3. Surface-pipeline models

TransX 0.00

24.00

48.00

72.00

96.00

Fig. 2. Reservoir description; heterogeneity and well placement.

2.2. Well vertical-flow models The vertical well flow model is integrated into the reservoir simulator by introducing the Tubing Head Pressure table (THP table). The THP table for each well is generated using the PROSPER® simulator. The data range and underlying model that are used to generate the THP table are provided in Table 3. The pressure-drop calculation for the production well in the gas-condensate reservoir is calculated using the Gray correlation (American Petroleum Institute Appendix B), Vertical Flow Correlation for Gas Well, while for oil reservoir is used Hagedorn and Brown correlation (Hagedorn and Brown, 1965). The application of the reservoir simulator to the well-reservoir system produces a tabulation of bottom-hole pressure versus surface rate, phase surface rate ratios, and tubing-head pressure. The data in the THP table reflects a particular PVT characterization, tubing size, length, roughness and geometric configuration. Bilinear interpolation is used to determine bottom-hole pressure for given values of rate, GOR (Gas Oil Ratio), GLR (Gas Liquid Ratio), water cut or WOR (Water Oil Ratio) and THP (Sensor Reference Manual, 2009). The producer rate constraint, the injector maximum bottom-hole pressure constraint and the plateau rate target are presented in Table 3.

The flow-line network connecting production wells to the firststage separator is represented by the surface-pipeline model. HYSYS® is used to calculate pressure losses in the pipelines. There are two types of pipelines, as shown in Fig. 3: one transports gas and the other liquid (oil and water). The pressure drop is solved through back-calculation; however, enough information must be supplied to complete the material and energy balance calculations. The solution procedure assumes that the pressure at the input to the 1st stage separator is known, and that the total field production rate, composition and temperature at the inlet of the pipeline, are known. HYSYS then performs a backward calculation to find the inlet pressure. The inlet pressure at the gas pipe is calculated using the Weymouth equation (Ikoku, 1984) and the pressure drop at the liquid pipe (oil and water) is calculated by using the Beggs and Brill correlation (Beggs and Brill, 1973). Heat transfer in the ground is assumed to be steady state and the same material is assumed in all pipes. The gas pipeline is assumed to be isothermal and the liquid pipeline non-isothermal. The pipeline data is presented in Table 4. 2.4. Surface-process model The surface model is a steady-state thermodynamic model where input streams vary with time because these inputs are determined by the reservoir models. The surface process model is implemented in HYSYS and is represented into two main separation processes: liquid and gas separations as shown in Fig. 3. The liquid separation process consists of multi-stage separation processes. Separators 1 and 4 are three-phase separation processes which separates gas, oil and water.

164

S.D. Rahmawati et al. / Journal of Petroleum Science and Engineering 81 (2012) 161–170 Table 4 Surface pipeline properties.

Table 3 Well and field properties. Parameter

Rate (units) Min:(intervals): max OGR

WGR

THP Tubing inside diameter Well depth Maximum producer rate constraint Minimum producer THP constraint Maximum injector BHP constraint Maximum plateau rate target

Parameter

Reservoir Lean GC

Rich GC

3

3

Sm /D 2831.69:(20): 1.42E+06 Sm3/Sm3 2.8E − 05:(10): 3.4E −03 Sm3/Sm3

Sm /D 2831.69:(20): 1.42E +06 Sm3/Sm3 2.8E −05:(10): 3.4E −03 Sm3/Sm3

0 Bara 6.89:(10):244.76 In 0.11 m 7425 Sm3/D 5.4 E +05 Bara 68.95 Bara 275.8 Sm3/D 2.7 E+06

0 Bara 6.89:(10):172.36 In 0.11 m 7425 Sm3/D 5.4 E +05 Bara 68.95 Bara 275.8 Sm3/D 2.7 E +06

Unit

Oil 3

Sm /D 15.9:(20): 3974.68 GOR (Sm3/Sm3) 53.4:(10):1781.1 Water cut (Sm3/Sm3) 0:(10):1 Bara 6.89:(10):344.73 In 0.13 m 8425 Sm3/D 1920 Bara 68.95 Bara 310.3 Sm3/D 9600

Length Inner diameter Roughness

km m mm

Reservoir Lean GC

Rich GC

Oil

5 0.254 4.60E −02

10 0.254 4.60E −02

11.5 0.3048 4.60E −02

Separators 2 and 3 are two-phase separation processes which separates gas and liquid. In sequence, the pressure for each separator is 56.2 bara, 21.7 bara, 4.5 bara and 1.01 bara. Further, there is a second-step drying stage for each separator to extract more liquid from the separated gas stream. The final product from the liquid separation process is condensate. A water pump is installed to transfer water to the water disposal facility. The gas separation process consists of CO2 removal, H2O removal and continues to the NGL plant; each process is simplified by representing it by a splitter model. In real field separation processes, they are represented as distillation columns. The Dew Point Controller (DPC) unit is installed to produce high NGL recovery. There are six final products from the gas separation process facility. These are sales gas, fuel gas, reinjected gas to the lean gas-condensate reservoir, reinjected gas to the rich gas-condensate reservoir, reinjected gas to the oil reservoir and NGL. There are two products from the NGL plant, NGL vapor and NGL liquid. NGL vapor mainly consists of methane and ethane and will be reinjected to the oil reservoir while NGL liquid mainly consists of heavy components which will be sold as NGL. Energy that is used or reproduced from unit operations is calculated

Fig. 3. Hysys surface process facility model.

S.D. Rahmawati et al. / Journal of Petroleum Science and Engineering 81 (2012) 161–170

as power consumption. In this integrated benchmark model, it is assumed that rigs are always available; product costumer demand, prevailing production prices, and HSE policy are always constant.

165

a) NPV surface response for different formulations and PTS, ( t p ) values. 7.0E+9 6.8E+9

2.5. Thermodynamic model

6.6E+9 6.4E+9

NPV (USD)

Peng–Robinson 1979 (Robinson and Peng, 1978 and Robinson et al., 1979) was used as the EOS (Equation Of State) model in the HYSYS, PROSPER and SENSOR simulators. A common EOS with 9 components was used to describe all gas-condensate reservoir fluids, another 6 components EOS was used to describe the oil reservoir fluids, and then the EOS models were continued into a “composite” EOS for surface calculation.

Annual Discounting (Eq. 4)

6.2E+9 6.0E+9 5.8E+9 5.6E+9

Project Time Step Discounting (Eq. 2)

5.4E+9 5.2E+9

2.6. Economic model

5.0E+9

 
Rc ðt Þ ¼ qg ðt Þr g þ qc ðt Þrc þ qNGL ðt Þr NGL − qwi ðt Þ þ qwp ðt Þ rw −pðt Þr p ð1Þ −M CO2 ðt Þr CO2 Þ

J NPV ¼

N X n¼1

"

# Rc ðt Þ −OPEX Δt p ð1 þ dÞn·Δt p ðyearÞ

RcNPV ðt Þ ¼ Rc ðΔt NPV Þ

ð2Þ

ð3Þ

200

400

600

800

1000

1200

Project Time Step (Days)

b) NPVs for a PTS (

t p ) of 1095 days with different calculation methods. 6.0E+06

Revenue & NPV (USD/D)

The goal of the integrated model is to study field asset value as represented by an economic model. This model is based on Net Present Value (NPV), with introducing a discount factor. The operational expenses (OPEX) cover the surface and well operation and are defined by a fixed amount, one million USD per day for the base case. The field revenue is obtained from gas (0.21 USD/m 3 (6 USD/Mcf)), NGL and condensate sales (503 USD/m 3 (80 USD/bbl)). The daily cost is summed from the volume of water production and injection (18.4 USD/m 3 (2.93 USD/bbl)), CO2 removal (15.4 USD/MT) and power consumption (5 cents/kWh). A key question is whether the Project Time Step (PTS) (Δtp) – i.e. the frequency of model integration update – has a significant impact on the NPV calculation. Frequent integrated-model updating can be expensive, as shown in Fig. 12 (lower curve). NPV consists of a net revenue stream, Rc(t), and a conversion to present value based on an annual discount factor d (e.g. d = 0.1/year). Formally, NPV can be calculated with any time discretization for a given revenue function, Rc(t). Convention is to calculate NPV annually, based on the average revenue during the year – even if the revenue changes significantly during the year. Rahmawati et al. (2010) converted Rc(t) to NPV with a time discretization equal to the PTS, and not using annual average discounting (Eqs. (1–2)). This leads to a misleading dependency of NPV on project time step, as shown in Fig. 4(a). The NPV kept increasing as PTS was decreased. This misleading behavior was found in the NPV calculation basis, as shown in Fig. 4(b). The red line in the figure shows the revenue in USD per day for a PTS 1095 days (3 years). The gray line represents an NPV calculation for each PTS while the black line is the NPV calculation on an annual basis. A more consistent discounting uses the same time discretization for NPV for all cases, independent of PTS (Eqs. (3–4)). This is also shown in Fig. 4(a), where it becomes clear that the project value, as expressed by NPV, has little dependency on PTS (Δtp). Rc ðΔt NPV Þ is the average annual cash flow, ΔtNPV is always 365 days and N = 20. The total simulation is 20 years.

0

NPV(DT-NPV=PTS) NPV(DT-NPV=365 Days) Revenue

5.0E+06 4.0E+06

Eq. 1

Eq. 4

3.0E+06 2.0E+06

Eq. 2

1.0E+06 0.0E+00 0

1000

2000

3000

4000

5000

6000

7000

8000

Time (Days) Fig. 4. NPV surface response.

3. Model based on decision support tools The I-OPT model is run by linking all software applications using the Pipe-It® software for integration and optimization. Pipe-It integrates and schedules SENSOR, HYSYS and Economic calculation for a given project run. Data transfer from one application to another provides dynamic communications among the simulators. HYSYS simulates the pipeline and surface-facility models and returns the injection compositions and injection rates to the SENSOR reservoir simulator through Pipe-It. The production rates, water injection rate, power consumption and mass of CO2 removal are transferred to the economic model. The compositional translation from the reservoir to the surface facility EOS characterization is solved by aggregating all components from the gas-condensate reservoirs and the oil reservoir. The algorithm for simulating the integrated model is shown in Fig. 5. The project time step and simulation end time must be preselected. The simulation ending time is used to define when the field operation is stopped. In the initial run, each reservoir model is run for 1 day. The reservoir simulation outputs are transferred to the surface model. The static surface simulation is also run for 1 day to obtain the gas injection compositions, gas injection rates, sales gas rate, NGL rate and condensate rate. The cash flow, Rc(t), is also calculated for this initial run. The next integrated model run is then controlled by the PTS (Δtp) value. 4. Base case simulation and sensitivity analysis

J NPV ¼

N  X R n¼1

 cNPV ðt Þ −OPEX Δt NPV ð1 þ dÞn

ð4Þ

The base-case data are shown in Table 5. The total simulation time is 20 years and injection is active during the first 10 years. The simulation

S.D. Rahmawati et al. / Journal of Petroleum Science and Engineering 81 (2012) 161–170

1.0E+10

9.0E+06

Define: 1. Project time step ( t p). 2. Simulation time end (tend)

Revenue (USD)

8.0E+06

Initialization of all model components For each project time step ( t p)

7.5E+09

Maximum NPV

7.0E+06

5.0E+09

6.0E+06

2.5E+09

5.0E+06

0.0E+00

4.0E+06

-2.5E+09 -5.0E+09

3.0E+06

Solve Reservoir Model

Compute the input to Hysys

Solve Hysys model

Calculate the revenue

PTS=1095 Days PTS=365 Days PTS = 30 Days

2.0E+06 1.0E+06

0

1000

2000

Cumulative NPV (USD)

166

-7.5E+09

3000

4000

5000

6000

7000

-1.0E+10 8000

Time (Days) If t = tend

Update: Gas injection composition, gas injection rate and t = t + tp

Fig. 6. Revenues and NPVs for different PTS (Δtp) using base-case parameters in Table 5.

gas increased after the end of the injection scenario, but later showed a downturn. Figs. 6 and 7 show that selection of PTS influences revenue performance. Hence, silo-model dynamics deserves further investigation. Here, we used four items to characterize silo-model performance: (1) Comparison of reservoir output, (2) Comparison of surface input, (3) Comparison of surface output simulations, and (4) Comparison of injection composition and rates. From the study, it can be concluded that:

Calculate the NPV based on annual basis Fig. 5. Numerical method for the integrated model.

scenario starts with injection for 10 years (gas injection for the gascondensate reservoirs and WAG injection for the oil reservoir), followed by depletion of the gas-condensate reservoirs, and water injection for the oil reservoir. The base-case WAG scenario is based on Scenario 2 SPE 5 Comparative solution project, (Killough and Kossack, 1987). The gas injection target rate is 566,336 m3/D, the water injection target rate is 7154 m3/D and the change from water to gas injection and vice versa occurs every 91.25 days. There are two active constraints for the oil reservoir, a GOR (1781 Sm 3/m 3) and a watercut (0.83). A well will shut down if it reaches one of these constraints and will re-open one year later. It may be noted that the water supplied for the water injection comes from an external source; hence, it is not directly linked to the process facility. The annual NPV performance for the base case is presented in Fig. 6. This figure shows that for the base-case parameters, the field should be operated for 10 years, from an economic point of view. A varying PTS (Δtp) does not change (annual-based) NPV significantly. As can be expected however high frequency dynamics are captured only for small PTS values. This is quite pronounced when observing the revenues for the first 10 year period for Δtp = 30 days. Please note that in Fig. 6 the NPV calculations are made with Eq. (4) and ΔtNPV = 365 days. Fig. 7 shows the sales gas, NGL, and condensate performances for the base-case parameters. This figure shows that sales

Table 5 Base-case parameters. Variable

Value

Sales gas fraction Fuel gas fraction Gas-condensate reinjection fraction Lean reinjection fraction DPC temperature (C) Discount factor (%) Injection time (days) Project time step (days) Total simulation time (days)

0.4 0.3 0.6 0.5 − 30 10 3650 365 7300

(1) Comparison of reservoir output. The oil reservoir produces heavier components (C3 +) than the two gas-condensate reservoirs. Here, the rich gas-condensate reservoir produces more (C3 +) than the lean gas-condensate reservoir, as shown in Fig. 8. There are some reservoirs that stop producing before 20 years, especially the gas-condensate reservoirs. Premature abandonment is caused by the reservoir pressure being too low to lift the liquid to the surface. Using a PTS of 30 days clearly shows dynamic reservoir output behavior caused by the composition changes inside the reservoir due to production and injection. (2) Comparison of surface input. The output from each reservoir consists of different components are combined into one input to the surface simulation, as it is explained in Section 3. It can be concluded from Fig. 9 that during injection period, C1 and C2 are dominantly produced but then it decreases when the injection is stopped.

3.5E+06

10500

3.0E+06

9000

2.5E+06

7500

2.0E+06

6000

1.5E+06

4500

1.0E+06

3000

5.0E+05

1500

0.0E+00 0

1000

2000

3000

4000

5000

6000

7000

0 8000

Time (days) Sales Gas - PTS 30D NGL - PTS 365D

Sales Gas - PTS 365D Condensate - PTS 30D

NGL - PTS 30D Condensate - PTS 365D

Fig. 7. Field products in volume units for different PTS values.

NGL and Condensate (Sm3/D)

t = tend

Sales Gas (Sm3/D)

No

0.35

1.4

0.3

1.2 1

0.25 0.2

0.8

0.15

0.6

0.1

0.4

0.05

0.2

0

0

1000

2000

3000

4000

5000

6000

7000

0 8000

Moles fraction (C3+/C2-) for oil reservoir

Moles fraction (C3+/C2-) for gascondensate reservoir

S.D. Rahmawati et al. / Journal of Petroleum Science and Engineering 81 (2012) 161–170

Time (Days) PTS30D - Lean PTS30D - Rich PTS30D - Oil

PTS365D - Lean PTS365D - Rich PTS365D - Oil

PTS1095D - Lean PTS1095D - Rich PTS1095D - Oil

Moles Fraction (C3+/C2-) from 3 reservoirs

Fig. 8. Reservoir outputs for different PTS values.

1 Hysys-Input (Total) PTS=30D Hysys-Input (Total) PTS=365D Hysys-Input (Total) PTS=1095D

0.9 0.8

167

The base case simulator and scenario were analyzed by perturbing several key parameters, i.e., the key decision variables in this benchmark case, including: (1) The dew point temperature controlling NGL extraction, (2) The gas sales fraction (fraction sales gas of total produced gas, TEE1 top-right in Fig. 3), (3) The gas-condensate reinjection fraction (fraction of reinjected gas into gas condensate reservoirs, TEE3 top-right in Fig. 3) and (4) The lean reinjection fraction (fraction reinjected gas into lean reservoir, TEE4 top-right in Fig. 3). For the reservoir aspect, it is possible to optimize the WAG period and the amount of gas and water injection rates. All other decision variables were held constant during a simulation run. Fig. 10 (a) is an example of single parameter analysis for an optimization variable and Fig. 10 (b) shows surface parameter analysis when two parameters are changed simultaneously. The parameter sensitivity results are summarized in Table 7 and Fig. 11. Fig. 11 shows minimum and maximum NPV ranges for different parameter sensitivities. Sensitivity of sales gas fraction and DPC temperature gives the highest maximum NPV, while sensitivity of WAG cycle period gives the lowest maximum NPV. Fig. 10 (a) and (b) show the day at which a maximum NPV is reached. It can be concluded from the study that nonlinear effects are significant and thus local optima are present in this problem. 5. Optimization

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

1000

2000

3000

4000

5000

6000

7000

8000

Time (days) Fig. 9. Surface simulation inputs from three different reservoirs.

(3) Comparison of surface output simulation. The surface output simulation in volume units is presented in Fig. 7. The sales gas is connected with the gas injected into the oil reservoir. When the oil reservoir is injected with water, the sales gas has a high value; when the oil reservoir is injected with gas, the sales gas decreases. The water and gas alternating injection depends on the reservoir injectivity index and PTS. It is clear that the dynamics are less complex after the injection scenario is stopped at 3650 days. Hence, the dynamic response is well captured for PTS equal to 365 days after the first 10 year period as opposed to the injection period. The cumulative production and revenues for different PTS values are presented in Table 6. The differences in cumulative undiscounted revenues among the PTS values are within approximately ~3%. (4) Comparison of injection composition and rate. The compositions of the injection products from the surface-process have a significant influence on reservoir production, as shown in Fig. 8.

The sensitivity analysis above indicates that optimizing NPV is nontrivial because the problem is non-convex. Our optimization is conducted by implementing the Nelder and Mead (1965) reflection simplex algorithm, modified to handle constraints and variable bounds (Pipe-It reference manual, 2010). This algorithm is a direct search method for nonlinear problems (Lagarias et al., 1998) where derivatives are not available or reliable. (Nocedal and Wright, 2006) explained that this algorithm seeks to remove the vertex with the worst function value and replace it with another point with a better value. The new point is obtained by reflecting, expanding or contracting the simplex along the line joining the worst vertex with the centroid of the remaining vertices. If the algorithm cannot find a better point in this manner, the algorithm retains only the vertex with the best function evaluation, and the algorithm then shrinks the simplex by moving all other vertices toward this value. The convergence criteria used in the optimization method are presented in Eqs. (5) and (6), representing NPV and variable movement tolerances. The optimization will converge if either criterion is met.        −3  J NPV best −J NPV worst ≤5·10 J NPV best  jxbest −xj≤10

ð5Þ

   xupper−bound −xlower−bound 

−6 

ð6Þ

The method is applied for three different optimization scenarios to maximize the NPV. The decision variables for the first scenario are DPC temperature, sales gas fraction, gas-condensate reinjection fraction and lean gas-condensate reinjection fraction. The decision variables for the second scenario are sales gas fraction, DPC temperature, gas injection rate and water injection rate for WAG scenario and WAG period. The decision variables for the third scenario are

Table 6 Cumulative production and revenue for different PTS. PTS (days)

CO2 removal (MT)

Sales gas (Sm3)

NGL (Sm3)

Condensate (Sm3)

Power consumption (KW)

Water injection & production (Sm3)

Cummulative revenue (USD)

30 365 1095

1.11E +07 1.18E +07 1.11E +07

1.11E +10 1.05E +10 1.16E +10

5.84E +06 5.52E +06 5.90E +06

3.61E+07 3.59E+07 3.61E+07

1.27E +08 1.21E +08 1.28E +08

4.12E +07 4.12E +07 4.12E +07

2.23E+10 2.20E+10 2.25E+10

168

S.D. Rahmawati et al. / Journal of Petroleum Science and Engineering 81 (2012) 161–170

a) Single parameter analysis for the WAG cycle time.

b) Surface parameter analysis for DPC temperature

This figure indicates that the best scenario for the oil reservoir is implementing Simultaneous WAG (SWAG).

and sales gas fraction. The maximum NPV is obtained when fsg = 0.5 and TDPC = -55 oC.

6.10E+09 3650

Base Case

6.00E+09

3650 4015 3650

5.95E+09

3650 3650 3650

5.90E+09

4015

5.85E+09 5.80E+09

Base Case

3650

6.8E+09

3650

Maximum NPV (USD)

Maximum NPV (USD)

3650

6.05E+09

4015

0

50

100

150

200

250

300

4015 4015

350

6.6E+09 6.4E+09 6.2E+09 6E+09 5.8E+09 5.6E+09 5.4E+09

400

-55 -50 -45 -40 -35 -30 0.1 -25 -20 -15 -10

0.9 0.7 0.5 0.3

DPC Temperature (C)

Delta WAG Cycle (days) Fig. 10. Sensitivity parameter analysis for single and surface parameters.

combination of the first and second optimization scenarios. The decision variables were defined based on the results from sensitivity study. The first scenario focused on the surface facilities parameter optimization, while the second and third scenarios are the combination of surface facilities and reservoir parameters. These optimization models can be described as the following:

Scenario 2 

max f sg ;T DPC ;qgi ;qwi ;Δt WAG

J NPV

O

with the following constraints on the decision variables 0:1≤f sg ≤0:9; −55≤T DPC ≤30; 0≤qgiO ≤1:81E þ 06; 0≤qwi

Scenario 1 n

max

≤8744:30; 30≤Δt WAG ≤365

o J NPV

f sg ;f Rgc ;f R ;T DPC

Scenario 3

L

with the following constraints on the decision variables 0:1≤f sg ≤0:9; 0:1≤f Rgc ≤0:9; 0:1≤f RL ≤0:9; −55≤T DPC ≤30

o J NPV

max

n

f sg ;f Rgc ;f R ;T DPC ;qgi ;qwi ;Δt WAG L

O

with the following constraints on the decision variables Table 7 Sensitivity parameter results. The optimal value and corresponding NPV are shown for each of the single parameter and two parameter sensitivity analyses, respectively. Case

DPC temperature Sales gas fraction Gas-condensate (GC) reinjection fraction Lean GC reinjection fraction WAG cycle period qgi and qwi for oil reservoir

GC inj fract and lean GC inj frac Sales gas frac and DPC temp

Optimal value

NPV Base case (USD) value

− 55 C

6.43E +09 6.17E +09 6.28E +09

0.3 0.1

0.6

30 days qgi = 8.5E + 05 m3/D and qwi = 5564.5 m3/D fRgc = 0.1 and fRL = 0.6

− 30 C 0.4 0.6

6.23E +09

0.5

6.09E +09 6.52E +09

91.25 days

6.43E +09

qgi = 5.7E + 05 m3/D and qwi = 7154.4 m3/D fRgc = 0.6 and fRL = 0.5

6.63E +09

fsg = 0.4 and Temp DPC = − 30 C

NPV – Increment (%) base case (USD) 6.03E +09 6.03E +09 6.03E +09

6.22

0:1≤f sg ≤0:9;

0:1≤f Rgc ≤0:9; 0:1≤f RL ≤0:9;

−55≤T DPC ≤30;

0≤qgiO ≤1:81E þ 06; 0≤qwi ≤8744:30;

30≤Δt WAG ≤365 The base-case parameters in Table 5 were used as the initial value for the optimization. The optimization results for scenario 1 are: fsg = 0.40, fRgc = 0.62, fRL = 0.69, and TDPC = − 55 oC. The optimization results for scenario 2 are: fsg = 0.1, TDPC = − 52.7 oC, ΔtWAG = 35.3days,

2.23 3.92

WAG Cycle Period Sales Gas Fraction

6.03E +09

3.21

6.03E +09 6.03E +09

0.90

Lean GC Reinjection Fraction Gas-Condensate Reinjection Fraction

7.56

GC Injection Fraction and Lean Injection Fraction

DPC Temperature Gas and Water Injection for oil reservoir

6.03E +09

6.22

6.03E +09

9.08

Sales Gas Fraction and DPC Temperature

4.00E+09 fsg = 0.5 and Temp DPC = − 55 C

5.00E+09

6.00E+09

NPV (USD) Fig. 11. Tornado chart of the parameter sensitivity in Table 7.

7.00E+09

S.D. Rahmawati et al. / Journal of Petroleum Science and Engineering 81 (2012) 161–170

6. Conclusions

Nomenclature d Discount factor N Total project time step MCO2 Mass of CO2 removal, MT/day p Power consumption, KW/day fsg Sales gas fraction fR Reinjected gas fraction, fR = 1 − fsg

Cumulative sales gas (m3) Cumulative NGL (m3) Cumulative condensate (m3) NPV (USD) Number of iterations CPU run time (hour) Increment from the base case

Base case

1.05E +10 5.52E +06 3.59E +07 6.03E +09 1 0.09 –

7.0E+9

30 25

Scenario 1 Optimum

6.5E+9

20

Base Case

6.0E+9

15

5.5E+9

10

5.0E+9

5

4.5E+9 0

200

400

600

800

1000

1200

Maximum NPV (USD)

7.5E+9

4.0E+9 1400

Project time step (days) Fig. 12. NPV for different project time steps.

fRgc fRO fRL fRR qgiO qwi qwp rc rNGL rp rw rCO2 TDPC Δtp ΔtNPV ΔtWAG

Gas-condensate reinjection fraction Oil reinjection fraction, fRO = 1 − fRgc Lean gas-condensate reinjection fraction Rich gas-condensate reinjection fraction, fRR = 1 − fRL Gas injection rate for oil reservoir, m 3/day Water injection rate for oil reservoir, m 3/day Surface water production rate, m 3/day Condensate price, USD/m 3 NGL price, USD/m 3 Power cost, USD/kWh Water injection and production cost, USD/m 3 CO2 removal cost, USD/MT Dew Point Controller temperature, °C Project time step, days Project time step equal to 365 days WAG cycle time, days

Acknowledgements We acknowledge our collegues Mohammad Faizul Hoda from PERA and Aleksander Juell from the Norwegian University of Science and Technology (NTNU). We are grateful for the assistance from Prof. Michael Golan at NTNU, Robert Hubbard from John M. Campbell & Co and Ken Starling from Starling & Associates, for helping with the surface process development. We also acknowledge the financial support by the Center for Integrated Operations in the Petroleum Industry at NTNU. References

Table 8 Comparison of base-case and optimization results. Parameter

Scenario 2 Optimum

35

0

The presented I-OPT model is suitable for assessing the potential of integrated optimization since the upstream and downstream parts of the model are tightly coupled. All aspects of the model are realistic and well suited for both life-cycle analysis and shorter time-frame studies. Detailed documentation is made available so alternative software platforms with the necessary functionality may be used to study the same multi-field, integrated asset system. The base case run time for the presented implementation on a standard laptop computer is ~ 6 min. Optimization has a clear potential because the multi-variable scenarios considered in this paper reveal an NPV increase of 10%–25% compared to the base case gas injection scenario. Fig. 12 shows that the magnitude of total NPV error is more-orless constant for a given Δtp(t), in both base-case and the optimization cases in Scenarios 1 and 2. We therefore concluded that the surface of maximum NPV is rather insensitive to Δtp(t), and thus compromised using a Δtp(t) of 1 year for the optimizations. Once an optimal case is located, the I-OPT project is rerun with a smaller project time step (e.g. 1 month) to obtain a more-accurate (“true”) value of the maximum NPV. Fig. 12 also shows that there was only a small gain to be made in terms of run time if the project time step is increased beyond 1 year. However, a shorter project time step increases the computational cost substantially. The NPV is shown for the varying project time steps. The base case NPV curve is equal to Fig. 4(a). One might argue for the selection of different project time steps depending on the run time and hardware resources available. Juell et al. (2010) improved the NPV result for a given project time step by introducing intermediate “division” project time steps whereby reservoir results were fed to the (fast and approximate) proxy surface process model without feedback. This approach was not used in our benchmark because the surface process is a real and more complex model, and therefore CPU time was much higher.

8.0E+9

40

CPU time (relative to the base case time)

qgiO = 1.74E + 06m3/D, and qwi = 3332.06 m3/D. The optimization results for scenario 3 are: fsg = 0.1, fRgc = 0.47, fRL = 0.59, TDPC = − 50.3 oC, ΔtWAG = 41.2 days, qgiO = 1.79E + 06 m3/D, and qwi = 3315.32 m 3/D. The comparison between the base case and the optimization results is presented in Table 8 and shows the potential of optimization since NPV has increased ~9%, ~23%, and ~25%, respectively. Scenario 1 required 29 iterations to converge to the optimum solution, scenario 2 took 73 iterations, and scenario 3 required 157 iterations. The base case CPU run time for optimization scenario 1 was ~2 h, it was ~7 h for Scenario 2, and it was ~12 h for Scenario 3.

169

Optimization Scenario 1

Scenario 2

Scenario 3

1.12E +10 6.75E +06 3.58E +07 6.61E +09 29 2.06 ~ 9%

1.11E +10 7.45E +06 4.66E +07 7.82E +09 73 7.02 ~ 23%

1.24E +10 7.58E +06 4.68E +07 8.00E +09 157 12 ~ 25%

American Petroleum Institute Appendix B, 1978. API 14 B. Vertical Flow Correlations – Gas Wells. p. 38. Bailey, W.J., Couet, B., Wilkinson, D., 2005. Framework for Field Optimization to Maximize Asset Value. SPE 87026-PA. SPE Reservoir Eng. J. 8 (1), 7–21. Beggs, D.H., Brill, J.P., 1973. A Study of Two-Phase Flow in Inclined Pipes. J. Pet. Technol. 607. doi:10.2118/4007-PA May. Cullick, A.S., Heath, D., Narayanan, K., April, J., Kelly, J., 2003. Optimizing Multiple-Field Scheduling and Production Strategy with Reduced Risk. Paper SPE 84239 presented at the 2003 SPE Annual Technical Conference and Exhibition, 5-8 October 2003, Denver, Colorado, U.S.A. doi:10.2523/84239-MS. Foss, B.A., Halvorsen, I.J., 2009. Dynamic optimization of the LNG value chain. Proceedings of the 1st annual gas processing symposium Doha, Qatar. Hagedorn, A.R., Brown, K.E., 1965. Experimental study of pressure gradients occurring during continuous two-phase flow in small-diameter vertical conduits. J. Pet. Technol. 17 (4), 475–484. Hepguler, G., Barua, S., 1997. Integration of A Field Surface and Production Network With A Reservoir Simulator. SPE 38937-PA. SPE Comput. Appl. J.. doi:10.2118/38937-PA Hepguler, G., Dutta-Roy, K., Bard, W.A., 1997. Applications of A Field Surface & Production Network Simulator Integrated With A Reservoir Simulator. SPE 38007-MS presented at SPE Reservoir Simulation Symposium, 8-11 June 1997, Dallas, Texas, U.S.A. doi:10.2523/38007-MS. Ikoku, C.U., 1984. Natural Gas Reservoir Engineering. Wiley, New York.

170

S.D. Rahmawati et al. / Journal of Petroleum Science and Engineering 81 (2012) 161–170

Juell, A., Whitson, C.H., and Hoda, M.F. 2009. Model Based Integration and Optimization Gas Cycling Benchmark. Paper SPE 121252. SPE Journal, September 2010. doi:10. 2118/121252-PA. Killough, J.E., Kossack, C.A., 1987. Paper SPE 16000. Fifth SPE Comparative Solution Project: Evaluation of Miscible flood Simulators. Presented at the Ninth SPE Symposium on Reservoir Simulation, 1-4 February 1987, San Antonio, Texas, U.S.A. Lagarias, J.C., Reeds, J.A., Wright, M.H., Wright, P.E., 1998. Convergence properties of the Nelder–Mead simplex method in low dimensions. SIAM J. Optim. 9 (1), 112–147. Lasschuit, W., Thijssen, N., 2004. Supporting supply chain planning and scheduling decisions in the oil and chemical industry. Comput. Chem. Eng. J. 28, 863–870. doi:10.1016/j.compchemeng.2003.09.026. Nocedal, J., Wright, S.J., 2006. Numerical Optimization Second Edition. Springer. Pipe-It Reference Manual, 2010. Zick, A. Reflection Solver. Rahmawati, S.D., Whitson, C.H., Foss, B.A., Kuntadi, A., 2010. Multi-Field Asset Integrated Optimization Benchmark. Paper SPE 130768 presented at the SPE EUROPEC/EAGE Conference, 14-17 June 2010, Barcelona, Spain.

Robinson, D.B., Peng, D.Y., 1978. The Characterization of the Heptanes and Heavier Fractions. Research Report 28. Gas Producers Assn., Tulsa, Oklahoma. Robinson, D.B., Peng, D.Y., Ng, H.Y., 1979. Capabilities of the Peng–Robinson Programs, Part 2: Three-phase and hydrate calculations. Hydrocarbon Proc. 58, 269. Sensor Reference Manual. Coats Engineering , Inc. 2009. Tomasgard, A., Rømo, F., Fodstad, M., Midthun, K., 2007. Optimization Models for the Natural Gas Value Chain. Geometric Modeling, Numerical Simulation, and Optimization. Springer, pp. 521–558. Trick, M.D., 1998. A Different Approach to Coupling A Reservoir Simulator With A Surface Facilities Models. Paper SPE 40001-MS presented at SPE Gas Technology Symposium, 15-18 March 1998, Calgary, Alberta, Canada. Yang, D., Zhang, Qi, Gu, Y., 2003. Integrated optimization and control of the productioninjection operation systems for hydrocarbon reservoirs. J. Pet. Sci. Eng. 37 (2003), 69–81.