- Email: [email protected]
Strategic optimization of the oil sands development with SAGD: Drainage area arrangement and development planning
Author’s Accepted Manuscript Strategic Optimization of the Oil Sands Development with SAGD: Drainage Area Arrangement and Development Planning Hossein Shahandeh, Shahed Rahim, Zukui Li www.elsevier.com/locate/petrol
PII: DOI: Reference:
S0920-4105(15)30187-X http://dx.doi.org/10.1016/j.petrol.2015.11.023 PETROL3256
To appear in: Journal of Petroleum Science and Engineering Received date: 24 June 2015 Revised date: 17 October 2015 Accepted date: 25 November 2015 Cite this article as: Hossein Shahandeh, Shahed Rahim and Zukui Li, Strategic Optimization of the Oil Sands Development with SAGD: Drainage Area Arrangement and Development Planning, Journal of Petroleum Science and Engineering, http://dx.doi.org/10.1016/j.petrol.2015.11.023 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Strategic Optimization of the Oil Sands Development with SAGD: Drainage Area Arrangement and Development Planning Hossein Shahandeh, Shahed Rahim, Zukui Li
Department of Chemical and Materials Engineering, University of Alberta, Edmonton, AB, T6G2V4 Canada
Abstract The majority of the oil sand deposits in Alberta Canada can only be extracted using in situ methods, mainly the Steam Assisted Gravity Drainage (SAGD). SAGD project consists of Central Processing Facility (CPF) and Surface Pad (SP) situated surface and Drainage Area (DA) consisting of multiple injector and producer well pairs situated subsurface. DA placement is crucial in ensuring that maximum amount of bitumen is extracted. In this work, an optimization framework is developed to plan the development of DAs. The first step aims to obtain a compact DAs arrangement by maximizing the amount of extractable bitumen. In the second step, a Mixed Integer Linear Optimization (MILP) model is developed to arrange the multiperiod development plan of the DAs. The proposed method is applied to a case study with multiple DAs for arrangement and the results demonstrate that the method can effectively generate a good DA layout and an economically optimal development plan that maximizes the Net Present Value (NPV).
Keywords Oil sands, SAGD, reservoir development planning, mixed integer linear optimization
Author to whom correspondence should be addressed. Email: [email protected]; Tel: 1-780-492-1107; Fax: 1-780-492-2881 1
1. Introduction Alberta’s oil sands are the third largest proven global reserves of oil with current estimates of 26.6 billion cubic meters of crude oil (Teare et al., 2013). Most of the oil sands in Alberta are discovered in three geological areas: Athabasca, Peace River and Cold Lake. Among them, Athabasca is the largest in terms of area. About 20% of the bitumen reserves are extractable using surface mining methods. The rest 80% are more than 70 meters below the ground and are extractable using in situ methods. Prominent in situ methods currently used in the oil sands industry are Cyclic Steam Simulation (CSS) and Steam Assisted Gravity Drainage (SAGD), whereas SAGD accounts for the most widely used in situ methods in Alberta. Nowadays, more than 16 SAGD projects are under operation in Alberta and many more are under development. Some examples of the prominent SAGD operations with significant bitumen production are: MacKay River and Firebag operations of Suncor Energy, Foster Creek and Christina Lake operations of Cenovus Energy and Jackfish operations of Devon energy (Teare et al., 2013).
SAGD operations use super-heated steam to decrease the viscosity of the bitumen in the reservoir. As a result, the bitumen can easily flow and be pumped up to the surface. As shown in Figure 1, SAGD wells consist of two horizontal wells which are in parallel. Steam is injected to the upper injector well and the melted bitumen flows to the lower producer well by gravity (Butler, 1991). The Bitumen from the producer well is pumped to the surface and transported to Central Processing Facility (CPF). Except CPF, the major SAGD facilities include Surface Pad (SP) and Drainage Area (DA). SP is the surface facilities from which multiple producer and injector well pairs are drilled. DA includes a set of parallel wells attached to the SP. CPF is a surface facility which processes the oil-water emulsion, recycles the water and produces the super-heated steam required by the injector wells. Steam is supplied from the CPF to each of the SP through pipelines. SAGD well pairs are drilled in areas with high volume of recoverable bitumen, which depends on the porosity and the oil saturation of the reservoir.
Current practices of commercial placement of SAGD wells largely rely on reservoir characterization and engineering judgement. Optimization methods for the placement of SAGDs’ DA and SP are relatively new in academia. Kumar investigated a space packing optimization method for generating a compact and non-overlapping arrangement of DAs (Kumar, 2011). Geometric transformations of global rotation, global translation and column translation were 2
used by the space packing algorithm to maximize the recoverable bitumen. Restrictions such as non-placement of SP over surface restrictions or non-placement of DA in thief zones were considered by incorporating a penalty function in the objective function. Manchuk and Deutsch used an adaptive grid search algorithm to determine the optimal arrangement of DAs (Manchuk and Deutsch, 2013). The optimization algorithm determined the positions and orientations of the SPs and DAs over a reservoir area to economically maximize the recoverable bitumen.
During the lifecycle of an oilfield project, development is the most critical phase due to intensive investments required for capital and operating costs. Determining whether it is economic to develop an oilfield or not and finding the best schedule for field development are the main concerns in this field of study. According to the timeframe, optimal decision making can be categorized as operational, tactical and strategic ones. On the other hand, exploration, development, production and transmission represent the major decision scopes (Shakhsi-Niaei et al., 2014).
Using simplifying assumptions such as linear performance of reservoirs versus time, optimization of the oilfield investment and operation planning was first tackled with Linear Programming (LP) (Aronofsky and Williams, 1962). Simultaneous optimization of locationallocation of wells, facility operation scheduling, and production rates during different time periods was then modeled as a Mixed Integer Linear Optimization (MILP) problem, but still based on linear performance of reservoirs (Frair, 1973). There have been many contributions in this field of study then after, and they can be mainly categorized into offshore and onshore oilfields’ development. Emphasizing on providing a detailed and accurate geological model of reservoirs, Behrenbruch outlined a systematic procedure for feasibility study of developing offshore petroleum fields (Behrenbruch, 1993). Tarhan et al. developed a multistage stochastic programming model for the planning of offshore oil/gas field infrastructure under uncertainty (Tarhan et al., 2009). Recently, Gupta and Grossmann proposed Mixed Integer Nonlinear Programming (MINLP) formulations for offshore oilfield development planning (Gupta and Grossmann, 2012; Gupta and Grossmann, 2014). In the former work, facility installation, well drilling and production decisions considering oil, water and gas flows profiles were simultaneously optimized when the problem was reduced to a MILP (Gupta and Grossmann, 2012). In the later work, a general multistage stochastic programming was proposed including
3
endogenous uncertainties (either fiscal rules or uncertainty in the field parameters) (Gupta and Grossmann, 2014).
Ierapetritou et al. proposed a novel decomposition-based approach for field development problem based on field decomposition and a quality cut-off criterion based iterative scheme (Ierapetritou et al., 1999). The optimization problem could be solved to optimality for each of the decomposed subproblem when it was formulated as a MILP problem. Kosmidis et al. presented a novel MINLP model considering nonlinear reservoir behaviour, multiphase flow in wells and constraints from the surface facilities (Kosmidis et al., 2005), by inclusion of piecewise linear approximations of each well model and using an outer approximation based algorithm. Tavallali et al. tried to apply a more realistic and detailed reservoir model and consequently efficient solution algorithm in two consecutive studies (Tavallali et al., 2014; Tavallali et al., 2013). The major contribution of the first work was the consideration of rigorous subsurface flow dynamics (Tavallali et al., 2013). They also introduced a unified model integrating subsurface, wells, and surface levels of an upstream production project. In the second work, a model was developed for location-allocation problem of well-drilling and infrastructure installation when the oilfield has more than one reservoirs with a shared surface processing network (Tavallali et al., 2014). Both the interactions within the surface network facilities and the dynamic states of the subsurface environment were incorporated into the model.
As the conventional oil resources become unavailable or uneconomic for extraction, bitumen has become more and more favorable. As an in situ recovery method, SAGD can be an excellent case study in the scope of strategic development optimization since its facility and recovery mechanism are very different from the conventional methods. Continuously injecting steam, having rising and depletion stages during bitumen recovery, and distinctive geological characteristics of each DA make this process a unique one. With this motivation, an optimization framework for determining the optimal placement and development order of SAGDs’ DA and SP facilities is proposed in this study. The objective of the DA arrangement optimization is to determine DA and SP locations which would provide maximum accessibility of the available bitumen resources. The objective of the development plan optimization is to maximize the Net Present Value (NPV). The optimal SAGDs’ DA and SP arrangement and development plan are obtained through a two-stage optimization framework. In the proposed method, the derivative free optimization solver NOMAD is used in Matlab to perform DA configuration optimization, 4
the mixed integer linear optimization solver CPLEX is used in GAMS to perform the development planning optimization.
The rest of the paper is organized as follows. The problem statement and the workflow for the proposed framework are explained in section 2. The proposed DA arrangement optimization method is presented in Section 3 followed with a case study. Section 4 presents the optimization model and results for reservoir development planning. The paper is concluded in section 5.
2. Problem statement
There are two objectives related to the strategic planning of the SAGD development. The first objective is to obtain a compact arrangement of DAs within the reservoir lease area to achieve maximum bitumen recovery. The second objective is to generate a plan for developing the DAs (i.e., determine the starting and ending time within the whole lifetime of the project) in order to reach the optimal NPV. A two-stage optimization method is proposed to achieve the two objectives. An overview of the two-stage optimization method is given in Figure 2 and it is summarized as follows.
The input to the first stage optimization includes the reservoir grid, geological data of the reservoir, any surface restrictions on the reservoir grid and the fixed dimensions of each DA and SP in the compact arrangement. A Matlab program is developed to calculate the total available bitumen with a given DA arrangement. A derivative free optimization method is used for the first stage optimization. The decision variables include the orientation and position of the arrangement of the DAs. The orientation is represented by the angle of the compact DA arrangement from the horizontal x direction. The position is represented by the translation of the compact DA arrangement along both the x and y direction. In the optimization model, a compact and non-overlapping arrangement of DA is considered so as to maximize the accessibility of the available resources. The DAs and SPs are in rectangle shape and two possible locations of SP are considered for each DA. The DA and SP dimensions are considered to be fixed in this study and as a result, fixed numbers of parallel injector/producer wells per DA are considered. The complexity of the optimization model is significantly reduced by fixing those parameters. The optimizer runs till the stop criteria of the optimizer (e.g., maximum number of iterations) is 5
satisfied. The results of the optimizer are the optimal orientation and location of the compact set of DAs. One or two candidate SP locations will be available for each DA in the compact arrangement depending on the surface restrictions and lease boundary.
The output of the first stage optimization provides input to the second stage optimization for generating the DA development plan. Oil production flowrates of all DAs are calculated in Matlab based on rigorous predictive model during DAs' lifetime, and they are used in the second stage optimization. A MILP optimization model is proposed for the second stage of the DA development planning. The MILP optimizer aims to identify the optimal DA development order that maximizes the NPV. This enables an optimal selection of DA for development in different time period. Different scenarios are also considered for the capacity of CPF in the second stage of optimization in order to identify the best capacity and capital investment for the CPF.
3. Drainage area arrangement optimization
Bitumen deposits are found in porous mediums underground between impermeable rocks. The top surface of the bitumen deposit below the impermeable section of rocks is defined as the Top Continuous Bitumen (TCB) and the bottom surface of the bitumen deposit above the impermeable section of rock is defined as the Bottom Continuous Bitumen (BCB). The BCB and TCB of a bitumen deposit are determined from data obtained from well logs or exploratory drills. Net Continuous Bitumen (NCB) is the total thickness between the TCB and BCB that meets certain minimum criteria for reservoir quality (Kumar, 2011) and it is calculated as the number of net cells which are above a threshold porosity and permeability value. For each DA, the available bitumen is the quantity of bitumen available for production along all the wells within the DA. In this study, the NCB value is used to quantify the available bitumen. Figure 3 illustrates how NCB is calculated along the elevation of the reservoir, where the black portions of the reservoir have lower porosity and permeability value and therefore they are not included in the determination of NCB. For the ith DA, the available bitumen is calculated by the following equation
Riavailable
NCB
cDAi
c
arc So
(1)
6
where Riavailable is the available bitumen, So is the oil saturation, NCBc is the net continuous bitumen for cell c. arc is the fraction of area of cell c inside the ith DA and its calculation is illustrated in Figure 4. Constant oil saturation So is assumed for the whole reservoir in this work. As a result, the total NCB value calculated as
cDAi
NCBc arc is a good indication of the
amount of bitumen available for extraction by the SAGD wells. Since the oil saturation So appears as a constant coefficient in (1), its value does not affect the DA configuration. To optimize the DA configuration and maximize the total available bitumen, a Matlab program was developed which calculated the total available bitumen as a function of the orientation and position of the compact set of DAs. The input of the program includes:
Grid size with the lease area boundaries and surface restriction on the grid,
Geological data for each cell in the grid,
DA dimensions, SP dimensions, DA-SP spacing,
Orientation (θ) and positions (∆x, ∆y) of the compact arrangement of DAs.
For each DA arrangement plan defined by (θ, ∆x, ∆y), the total available bitumen is evaluated based on the following principles:
A DA is considered valid only if its area is 80% or more inside the lease boundary. For a rectangle shape DA, we assume it is worth to be developed only when the majority of its area is within the lease boundary. Otherwise, the major part of the DA will be outside the lease boundary such that its development with horizontal well is not practical.
A DA is considered valid only if its total available bitumen is above a certain threshold value. This assumption is made to avoid placing a DA in the region that has very low available bitumen and is not economical to be developed. This threshold is determined using the following method. Initially, a DA configuration with zero angle and zero translation is evaluated to find the distribution of the available bitumen for all the DAs. Based on the distribution, the threshold value for total available bitumen is set as the number such that 10% of the values lie below it.
A DA is considered valid only if it has at least one possible location for placing a SP. For a feasible SP location, the SP had to be placed inside the lease boundary and could not be placed on top of any surface restrictions (e.g., rivers, lakes, roads or conserved forest areas).
The output of the Matlab program is the total available bitumen from the given DA arrangement. 7
The first stage in the optimization framework generates a compact arrangement of DA. In this step, the objective function for the compact DA arrangement optimization problem is designed as maximizing the total available bitumen, which can be represented as
max Riavailable ( , x, y )
(2)
i
where Riavailable is the available bitumen for the ith DA. However, there is no explicit equation for the above objective function since it is evaluated through a Matlab program. Derivative free optimization technique is used in this work which calls on the developed Matlab program for total available bitumen calculation. Specifically, the derivative free optimization solver NOMAD is used, which implements the Mesh Adaptive Direct Search (MADS) algorithm for constrained blackbox functions. The MADS algorithm is an extension of pattern search method for nonlinear constrained optimization problems (Audet et al., 2009). During this optimization, upper and lower bounds were applied for the decision variables: orientation, translation along x direction and translation along y direction (by assuming the left bottom corner of the reservoir grid is initially located at the origin with x direction parallel to the horizontal axis), as given by Equations 3–5 respectively.
(3)
x x x
(4)
y y y
(5)
where is radian, x and y are twice the length and width of the reservoir, respectively.
Case Study
A case study with a three dimensional reservoir grid is used to demonstrate the practicality of the optimization method used for SAGD DA placement. A relatively large reservoir model with a grid size of 55 × 62 × 10 cells with cell dimensions of 100m × 100m × 5m is considered. The geological data considered for the case study are the porosity and permeability of the reservoir. Porosity and permeability were generated for every cell of the grid using a built-in function GaussianField in Matlab Reservoir Simulation Toolbox (Lie et al., 2012). The porosity values were in the range of 0.1 to 0.5. Permeability was calculated from the porosity data using Carmen-Kozeny relationship. Furthermore, in this case study an irregular shaped leased 8
boundary was used to reflect the common situation of oil sands projects. In addition, surface restrictions were created within the lease boundary to reflect surface restrictions for placing SP. The aerial view of the lease boundary and the surface restriction is given in Figure 5. In Figure 5, the leased area is represented as the yellow section, surface restriction is represented by the blue area and the white parts represent the area outside the lease boundary. Notice that the surface restriction area is virtual in this case study, and it has been divided into convex polygons for use in the Matlab code of NCB evaluation. Every SP is forced to have no interaction with any of those polygons. However, DA can be placed under those surface restrictions since it is underground.
The porosity and permeability distribution of the three dimensional grid is given in Figure 6 and 7 respectively. The dark blue regions in those figures represent the unleased area. Geological data for the unleased area is considered to be zero to simplify the recoverable bitumen calculation.
The geological data is used to calculate the NCB of the reservoir. In this case study, any cell with a porosity and permeability above a threshold value of 0.25 and 1.1×10-13 m2 respectively, were considered to be a net cell. A threshold value for the total bitumen is applied for each DA to be counted in the objective function (2). It is determined as 300 in this case study using the aforementioned method. The aerial view of the NCB of the reservoir grid is given in Figure 8.
In order to provide SAGD well arrangement plan, a list of parameters specific to all the DAs and SPs in the compact arrangement were fixed. The fixed dimensions of the DA and SP are given in Table 1. Based on the setting, each DA has six pairs of injector and producer wells. Two possible locations of SP with respect to a specific DA were considered with a fixed DA-SP spacing.
Overall placement plan The proposed optimization method was applied to the reservoir and it provides a compact arrangement of DAs. The results of orientation and translation along the x and y direction for the compact DA layout is given in Table 2. The total available bitumen in Table 2 is the sum of the available bitumen values of all the DAs in the compact arrangement.
9
The compact arrangement of all the DAs from the first stage optimization is given in Figure 9. The numbers on each of the DA quantifies the amount of available bitumen (NCB value) from that particular DA. The following observations can be made from Figure 9:
No DA is placed outside of the leased area
No DA is placed on areas which would result in a significant fraction of DA area to be outside the lease boundary
No SP is placed outside of the leased area
No SP is placed on any of the surface restrictions
It is important to note that due to surface restriction and lease boundaries, some DAs have only one possible location for the placement of SP.
Decomposition based placement The first stage optimization method was further studied by decomposing and breaking down the lease area into subareas. The DA configuration optimization was applied to each of those subareas to determine the optimal orientation and position for the compact set of DAs. Two lease area decomposition schemes are investigated here. The first scheme is shown in Figure 10, where the red dotted lines divide the lease area into 3 subareas. Table 3 lists the optimal orientation and position of the different DA arrangement for 3 subareas. The total available bitumen of each of the subareas and the sum of the total available bitumen from all the subareas are also reported in Table 3. The results show that no DA or SP is placed outside the lease boundary and no SP is placed on any of the surface restricted areas. Furthermore, a DA was not placed if the DA did not have any available position for placement of at least one SP. For the above decomposition scheme, the total available bitumen from the 3 divided areas (15,681) is lower than the total available bitumen from the undivided compact DA arrangement plan of Figure 9 (16,720).
The second decomposition scheme is shown in Figure 11. The DA arrangement optimization was then performed on the 3 decomposed sections and the total available bitumen from the 3 different sections are reported in Table 4. In this case, the total available bitumen from summation of the total available bitumen from all the different areas (16,947) is higher than the total available bitumen from the DA arrangement obtained from the undivided lease area (16,720). The above study shows that decomposing a large lease area into small subareas and then applying the compact DA arrangement optimization method separately to those subareas may result in higher total available bitumen. However, it requires testing different decomposition 10
scheme and finding the best way to divide the lease area. This can be deemed as another optimization problem for seeking the global optimal DA configuration. Consider the complexity of existing optimization problem, this study is not included in this work.
4. Reservoir development planning
The objective of the second stage optimization of the proposed framework is to determine the DA development plan based on the compact arrangement solution obtained from the first stage. In this section, the oil production rate profile is first calculated through simulation of nonlinear production models for each DA according to its geological specifications. Parameters used in the calculation are provided and explained. A MILP problem is then formulated for planning the DA development with the objective of maximizing the NPV. The DA layout shown in Figure 9 was used as the basis for the case study of this section.
4.1 Well production model
The crucial part of the second stage optimization is a model for oil production rate through SAGD process based on its geological characteristics. Some models have been developed for this purpose in which the SAGD production is divided into rising steam chamber stage and lateralspreading steam chamber stage (Akin, 2006; Guo et al., 2013; Vanegas et al., 2011), and the oil production rates in the two different stages are calculated as
akgSoh t 2 kg / m s qd 2 2 2 m s w bSoh 1/2
3/2
4 qr 3
9 42
2/3
(6)
2/3
So
1/3
kg 1/3 t m s
(7)
where qd and qr (m2 s-1) are oil production rate per meter of well length in depletion and rising stage respectively, k (m2) is the effective reservoir permeability, g (m s-2) is the acceleration due to gravity, α (m2 s-1) is the reservoir thermal diffusivity, ϕ is the porosity, ΔSo is the displaceable oil saturation, h (m) is the effective drainage height, m is a viscosity-temperature correlation parameter, νs (m2 s-1) is the oil kinematic viscosity at the steam temperature, t (s) is the time, w (m) is the half well spacing, γ is the shape factor of the rising steam chamber, and a, b, β are
11
empirical constants. Among the mentioned characteristics, ΔSo is the only one varying with time (Guo et al., 2013);
h So So 0.43 s kgt
0.4
(8)
where So is initial oil saturation of the reservoir. To predict the oil production rate, the depletion model and rising model are combined together with weighting factors (Guo et al., 2013):
qtoil yd qd yr qr
(9)
yd
1 1 10ud ( vd t )
(10)
yr
1 1 10ur ( vr t )
(11)
where qtoil is the total oil production rate, yd is the weighting factor of depletion model increasing from 0 to 1, yr is the weighting factor of rising model decreasing from 1 to 0, t (year) is the time, and ud, vd, ur, vr are constants. In the above equations, there are many nonlinear terms (especially for depletion stage) that will cause difficulty for optimization. Therefore, they were not directly used in the optimization model for oil production rate calculation. Instead, oil production rates were first calculated versus time for each DA, and the rates were then utilized as parameters to describe the oil production rate at each year. In this way, the optimization problem could be simplified and then solved with much less computational efforts.
Athabasca reservoir characteristics adapted from reference (Guo et al., 2013) were used to validate the presented model for oil production rate and also demonstrate trend of qd, qr, and qtoil versus time. Related parameters are provided in Table 5. Figure 12 shows the rising stage is the dominant stage at the initial stage and the depletion stage is the dominant one at the end, while they are both partially effective in the middle.
For the case study of this work, oil production rates of each individual DA can be predicted according to its reservoir characteristics. Average porosity, permeability, and initial oil saturation
12
are the parameters which can be defined for each DA based on the following equations (Rahim et al., 2015);
0.25 NCB DAcell
Av k
(12)
6 dp
(13)
1 3 2 Av 2 1 2
100 2.25 k Sw
(14)
2
(15)
So 1 Sw
(16)
where DAcell is the number of cells in each DA, Av (m-1) is the specific area of spherical shape grains, dp (m) is the average diameter of grains, τ is tortuosity, and Sw is water saturation. Constants are reported in Table 6. It should be highlighted here that each DA was assumed to have six pairs of injector and producer wells.
4.2 Planning optimization model
Based on the above well production model, the DA development planning problem can be defined next. In the model, indices i, t and k correspond to DA, time period within the lifetime of the whole project, and time period within the lifetime of a DA, respectively. The total number of DAs was 43 and the total lifetime of the project was set as 25 years. The total lifetime of each DA, DAiLife , was obtained from Matlab calculation and it varies from 6 to 10 years depending on its geological property. Note that 10 years is an upper bound of the life time for all the DAs based on the production rate calculation.
In the optimization model, the starting point ( tiStart ) and ending point ( tiEnd ) of each DA are the main decision variables. It is assumed that decisions can only be made at the beginning of each year and any operation of a DA can only be changed at that time. Starting point is defined by following equation: tiStart xi ,t t i
(17)
t
13
where xi,t is a binary variable which takes value 1 if ith DA starts its operation at time t, and value 0 otherwise.
To make sure that there is only one starting point for each DA, the next constraint is applied. If the summation of binary variables is equal to zero, it means that DA should not be operated at all during the lifetime of the project.
x
i ,t
1 i
(18)
t
The lifetime of each DA, tiLife , is a parameter varying from 6 to 10 based on its geological property. The ending point of a DA can be then determined as
tiEnd tiStart tiLife
x
i ,t
i
(19)
t
If being under operation of a DA is uneconomical, corresponding starting and ending times will be both equal to zero. There should be some mathematical logics between the oil production rate qioil,k of ith DA at kth year of its lifetime and the oil production rate of the same well in 25 years basis of the whole project f i oil ,t . In the other words, the lifetime of each DA can be shifted through the project horizon. This can be achieved by introducing a new binary variable yi ,k ,t . For each DA i, therefore, there exists a matrix of binary variables with k rows and t columns which can only have a value of one in a diagonal form. To determine the value of one in the first row of each matrix a big-M formulation can be applied and the oil production rate qioil,k would be equal to f i oil ,t if only t is equal to tiStart k 1 . The following equation was applied to reflect this logic: 1 yi ,k ,t M t k 1 tiStart 1 yi ,k ,t M t k 1 i, t, k 1
(20)
Equation 20 is a “big-M” constraint. It is active when yi,k,t is equal to 1, otherwise, it is redundant. Based on the inequality, M should be bigger than the largest possible value of
tiStart (t k 1) and tiStart (t k 1) , which would be 27 (=19 – 1 + 10 – 1) according to the variables’ domain: a well starting lifetime cannot be later than year 19th and its longest lifetime is 10 years. While any number that is bigger than 27 should be a valid option for M, we choose 28 as the value to keep the constraint tight.
14
For the rest of the rows, the remaining y binary variables can only be equal to one in a diagonal line,
yi ,k ,t yi ,k 1,t 1 i, t 1, k 1, k c1
(21)
As mentioned earlier, DAs have different lifetime varying from 6 to 10 years, but all the y matrices have 10 rows for the sake of mathematical consistency. As a results, 10th, 9th–10th, 8th– 10th, and 7th–10th rows of y matrices must be fixed to zero when the lifetime of corresponding DA is respectively 9, 8, 7, 6 years to avoid any incorrect mathematical logic. The c1 constant number in equation 21 is determined based on this concept. Moreover, there are upper and lower triangles in y matrices which also can set to zero due to diagonal form of binary variable equal to one. In this way, the number of binary variables can be reduced by more than 40%.
With the binary variables yi ,k ,t , the following equation can be used to shift any DA with different lifetime through the project horizon for calculating the oil production rate
yi,k ,t qioil,k k
fi oil ,t L DAwell
i, t
(22)
where L is the length of each well and DAwell is the number of producer wells under operation for each DA. Even though the left hand side of the last equation has a sigma sign, only one of the multiplication will be nonzero due to the way that yi,k,t has been defined. Note that qioil,k , which is calculated through equations 6–11, has the unit of (m3 m-1) and it only corresponds to the 3 flowrate of single pair of wells per its length. On the other hand, unit of f i oil ,t is (m ) and it shows
the whole flowrate of each DA including its six pairs of wells. The above three equations work properly if only the following one is also incorporated. Otherwise, all the yi,k,t would be equal to zero.
y
i ,k ,t
k
tiLife
i
(23)
t
The rest of the model includes economic model of the process which is based on NPV. It was assumed that the main utility consumption of oil production is related to steam. Required steam flowrate and its cost can be estimated as; fi ,steam fi oil t ,t SOR i , t
(24)
TFCt FC fi ,steam t t
(25)
i
15
where f i ,steam is required flowrate of steam for ith DA at tth time period, SOR is a constant number t used for Steam to Oil Ratio, TFC is Total Fuel Cost (M$ yr-1) of the tth year , and FC is Fuel Cost (M$ m-3).
Capital cost of process is composed of CCCPF (M$), initial investment for CPF which is mandatory, and CCDA (M$), related to investment for each DA installation. The value for CCCPF was estimated from the reference (Giacchetta et al., 2015) and a CPF capacity value of 3.5×106 m3 yr-1 was used as the base value, and all the 50%, 100%, 150%, 200%, 250%, and 300% of it were studied as possible scenarios. A proportional relation was implemented between the capital cost of CPF and its capacity. It was also assumed that each DA has 6 pairs of injector and producer wells and capital cost of each of them was 78 M$ (JacobsConsultancy, 2012).
TCCt CCCPF CCDA xi ,t
t 1
(26)
i
TCCt CCDA xi ,t
t 1
(27)
i
Equation 26 is only valid for the first year of project, and equation 27 is valid for the remaining ones. The only difference between equations 26 and 27 is that capital cost of CPF should be just taken into account at the first year. In the following equations, TOMC is Total Operating and Maintenance Cost (M$ yr-1), TAC is Total Annual Cost (M$ yr-1), TAR is Total Annual Revenue (M$ yr-1), OI is Operative Income (M$ yr-1), they can be calculated as (Giacchetta et al., 2015; Lawal, 2014);
TOMC
Discount rate TCCt t
(28)
SAGDlifetime
TACt TFCt TOMC t
(29)
TARt 1 Royalty OP qi ,t t
(30)
i
OI t TARt - TACt t
(31)
where OP is oil price (CAN$ m-3). If a company wants to extract oil sands in an area which is predefined as oil sands area by the government of Alberta, they have to pay a percent of gross revenues as royalty. This number has scaled value according to the price of West Texas Intermediate.
16
Tax can also be calculated using following equation (Giacchetta et al., 2015; Lawal, 2014);
Taxt TOI t Tax rate t
(32)
TOI t OI t Akt
(33)
Akt
t
xt CCCPF CCDA i ,Life SAGDlifetime i DAi
Akt CCDA i
xi ,t DAiLife
t 1
t 1
(34)
(35)
where TOI is Total Operating Income and Ak is amortization of investment.
NPV is calculated based on NI (Net Income) and DCF (Discount Cash Flow). They can be calculated as (Giacchetta et al., 2015; Lawal, 2014):
NI t OI t Taxt t
DCFt
(36)
NI t
1 Discount rate
t
t
(37)
NPV DCFt TCCt t SAGDlifetime t
(38)
t
Two additional constraints were imposed in this optimization problem. The first one was on the total amount of DAs installation’s investment at each year which should be equal or less than 300 M$.
TACt TCCt 300 t
(39)
Furthermore, total oil production at each year must be equal or less than the capacity of CPF. When a capacity is selected for a plant during design phase, processing of raw materials more than that fixed number would be impossible unless using the existing equipment more effectively or purchasing new equipment. While increasing the capacity has not been considered here as an option, the following constraint could be used as the upper limit of CPF capacity:
f
oil i ,t
3.5 106 c2 t
(40)
i
where c2 constant is a coefficient to consider different scenarios of 50%, 100%, 150%, 200%, 250% and 300% of the nominal capacity (3.5×106 m3/year) of CPF.
Finally, the optimization model can be expressed as 17
max NPV s.t. eqs(17 40)
(41)
The optimization problem (41) is a MILP problem. With all the required constants provided in Table 7, the model was solved in GAMS using the CPLEX solver. The optimization stopping criteria were set as one hour computational time and relative optimality gap of 1%. The solution for the base capacity of CPF (i.e., c2=1) was found in less than 804 seconds with a desktop computer (single core of Intel® i5-4590 @ 3.30 GHz, 8 GB RAM). The optimal NPV value at the end of the project lifetime was 12559.6 M$.
4.3 Result and discussion
As mentioned in previous section, 6 different scenarios for the capacity of CPF were taken into account to investigate the design characteristic of the project. Computational statistics for each case study are available at Table 8. To have a detailed study, the results of three capacity scenarios (50%, 150% and 200%) are illustrated through Figures 13–15. Each of them consists of a development plan (top left), a profile of the number of operating DAs (top right), NPV (bottom left) and total oil production (bottom right) over 25 years of project lifetime. Note that the segments in the development plan figure represent the operating years of each DA.
In this study, all the 43 DAs shown in Figure 9 were first sorted based on their NCB (bitumen) values in an ascending order. Therefore, the DA number 1 has the lowest NCB and the DA number 43 has the highest NCB. As the Figures 13-15 show, the DAs with higher NCBs are more favorable for being developed and under operation at the early stage of the project lifetime. When the capacity of CPF is small, many DAs with small NCBs are not even applied at all. On the other hand, when the CPF has a large capacity, all the DAs tend to be developed while the priority of operation is still on the DAs with high NCB values. If the DA’s investment constraint and CPF capacity constraint had not been imposed, all the DAs would have been developed and operated since the first year. The reason is that keeping the operation running for a longer time means higher amortization and operating cost. Activation of the investment limit constraint and CPF capacity constraint forces the project to take longer time of operation, and only three DAs could be installed during each year. Another interesting fact about these results is that the DA’s investment constraint was active at c2=50%, 150%, 200%, 18
while the constraint associated with the capacity of CPF was completely active at c2=50%, partially active at c2=150%, and completely inactive at c2=200%. Due to CPF investment in the first year, NPV starts to grow from a negative value. Afterwards, it becomes higher and higher versus operating year since the capital cost of DAs’ installation is not that much compared to the oil production revenue. In the last years, NPV does not change significantly since the production rate decreases more and more and there is not any investment for a new DA in the last five years.
The larger the CPF capacity, the higher number of DAs could be operated and also the higher total oil production could be achieved. Nevertheless, larger CPF requires higher capital cost investment which will lead to a lower NPV. As demonstrated in Figure 16, there exists an optimal point for the capacity of CPF which happens at 5.25 Mm3 yr-1 (150% of the base capacity).
Since the optimal scenario was found at this point, it would be beneficial to provide more information about it such as planning development of the case study on compact DA placement plan (Figure 17) and flowrate of each individual DA during project lifetime (Figure 18). Figure 17, which is based on the DA layout in Figure 9, shows different DAs and corresponding NCB (black number) and lifetime (yellow number) to them. For example, the DA at the top left corner has the NCB of 309 and is under operation from year 14th to 24th. It worth to mention that darker rectangular represent DAs with higher NCB and lighter ones represent DAs with lower NCB.
Oil production of each DA versus operating year is also provided in Figures 18. As it can be observed, DAs with better performance are applied for operation at beginning years. DAs’ oil production profiles are represented in Figure 18 and darker lines belong to DAs with larger ID (i.e., higher NCB values).
5. Conclusion
A novel optimization framework for strategic planning of in situ oil sands production has been proposed in this paper. The proposed method determines the compact arrangement of SAGD’s DAs and the optimal development order for those DAs. The first stage of the proposed 19
optimization method determines the optimal orientation and position of a compact set of all the DAs and SPs that can be packed within the lease area of the reservoir so that the total extractable bitumen is maximized. Surface restrictions and irregular shaped lease areas are taken into consideration in the placement of DAs and SPs. The SAGD arrangement optimization model produced a compact DA arrangement plan with high recoverable bitumen resource. The second stage of the proposed optimization method determines the optimal DA development plan for the lifetime of the project. This model generates a plan that maximizes the NPV under the investment and CPF capacity restrictions. Different scenarios were also studied for capacity of CPF to show the importance of this designing factor. The optimization method was applied to a case study with realistic reservoir size to determine a DA placement and development plan that maximizes the NPV from the SAGD project.
There are several ways to further improve the proposed method. First, the proposed DA configuration optimization problem is based on fixed dimensions of each DA. To achieve global optimal placement, non-uniform size of each DA should be modeled although it leads to a more complex optimization problem since the number of decision variables significantly increase. In addition, developing a optimization model for identifying the best possible decomposition of lease area is also a way to further improve the solution. Finally, in future work, the SAGDs’ DA arrangement optimization method will be studied by incorporating geological uncertainty (i.e., multiple geological realizations). Long term uncertainty in crude oil demand and price will also be considered to study the sensitivity and robustness of the reservoir development plan.
Acknowledgements The authors would like to acknowledge the financial support from Natural Sciences and Engineering Resource Council (NSERC) of Canada Discovery Grant Program.
20
Nomenclature Indices i
drainage area (DA)
c
cell
k
time period within the lifetime of a DA
t
time period within the lifetime of the whole project
Parameters a, b, β
empirical constants in oil production rate equation in depletion stage
Av
specific area of spherical shape grains (m-1)
CCCPF
initial investment for CPF (M$)
CCDA
investment for each DA installation (M$)
DAcell
number of cell in each DA
DAlife
total lifetime of each DA
DAwell
number of producer wells under operation for each DA
dp
average diameter of grains (m)
FC
Fuel Cost (M$ m-3)
g
acceleration due to gravity (m s-2)
h
effective drainage height (m)
k
effective reservoir permeability (m2)
L
length of each well (m)
m
viscosity-temperature correlation parameter
M
big-M constant
OP
Oil price (CAN$ m-3)
So
initial oil saturation of the reservoir
SOR
Steam to Oil Ratio
Sw
water saturation
ud, vd
empirical constants for weighting factor equation of depletion model
ur, vr
empirical constants for weighting factor equation of rising model
w
half well spacing (m)
α
reservoir thermal diffusivity (m2 s-1)
γ
shape factor of the rising steam chamber
νs
oil kinematic viscosity at the steam temperature (m2 s-1)
τ
tortuosity
porosity
Continuous variables Ak
amortization of investment 21
ar
fraction of area of cell c inside the ith DA
DCF
Discount Cash Flow
foil
oil production rate in 25 years basis of the whole project (m3 s-1)
fsteam
required flowrate of steam
NCB
Net Continuous Bitumen
NI
Net Income
OI
Operative Income (M$ yr-1)
qd
oil production rate per meter of well length in depletion stage (m2 s-1)
qoil
oil production rate in lifetime years basis of each DA (m2 s-1)
qr
oil production rate per meter of well length in rising stage (m2 s-1)
qtoil
total oil production rate (m2 s-1)
Ravailable
available bitumen for each DA
t
time (s)
TAC
Total Annual Cost (M$ yr-1)
TAR
Total Annual Revenue (M$ yr-1)
tend
ending point
TFC
Total Fuel Cost (M$ yr-1)
tlife
lifetime of each DA
TOI
Total Operating Income
TOMC
Total Operating and Maintenance Cost (M$ yr-1)
tstart
starting point
yd
weighting factor of depletion model
yr
weighting factor of rising model
ΔSo
displaceable oil saturation
∆x
horizontal position of the compact arrangement of all DAs
∆y
vertical position of the compact arrangement of all DAs
θ
orientation of the compact arrangement of all DAs
Binary variables xi,t
whether the operation of ith DA starts at time t or not
yi,k,t
whether the operation of ith DA starts at time k (DA lifetime) and time t (project lifetime) or not
22
References
Akin, S., 2006. Mathematical modeling of steam-assisted gravity drainage. Computers and Geosciences, 32(2): 240-246. Aronofsky, J.S. and Williams, A.C., 1962. The Use of Linear Programming and Mathematical Models in Underground Oil Production. Management Science, 8(4): 394-407. Audet, C., Digabel, S.L. and Tribes, C., 2009. NOMAD user guide, Les Cahiers du GERAD. Behrenbruch, P., 1993. Offshore oilfield development planning. JPT, Journal of Petroleum Technology, 45(8): 735-743. Butler, R.M., 1991. Thermal recovery of oil and bitumen. Prentice Hall Inc., Englewood Cliffs, NJ. Frair, L.C., 1973. Economic Optimization of Offshore Oilfield Development, University of Oklahoma, Tulsa, OK. Giacchetta, G., Leporini, M. and Marchetti, B., 2015. Economic and environmental analysis of a Steam Assisted Gravity Drainage (SAGD) facility for oil recovery from Canadian oil sands. Applied Energy, 142: 1-9. Guo, J., Zan, C., Ma, D.S. and Shi, L., 2013. Oil production rate predictions for steam assisted gravity drainage based on high-pressure experiments. Science China Technological Sciences, 56(2): 324-334. Gupta, V. and Grossmann, I.E., 2012. An efficient multiperiod MINLP model for optimal planning of offshore oil and gas field infrastructure. Industrial and Engineering Chemistry Research, 51(19): 6823-6840. Gupta, V. and Grossmann, I.E., 2014. Multistage stochastic programming approach for offshore oilfield infrastructure planning under production sharing agreements and endogenous uncertainties. Journal of Petroleum Science and Engineering, 124: 180-197. Ierapetritou, M.G., Floudas, C.A., Vasantharajan, S. and Cullick, A.S., 1999. Optimal location of vertical wells: Decomposition approach. AIChE Journal, 45(4): 844-859. JacobsConsultancy, 2012. Assessment of Innovative Applications of Electricity for Oil Sands Development Phase 1 Report. www.ptac.org/attachments/1425/download. Kosmidis, V.D., Perkins, J.D. and Pistikopoulos, E.N., 2005. A mixed integer optimization formulation for the well scheduling problem on petroleum fields. Computers & Chemical Engineering, 29(7): 1523-1541. Kumar, A., 2011. Optimal drainage area and surface pad positioning for SAGD development, University of Alberta, Canada. Lawal, K.A., 2014. Economics of steam-assisted gravity drainage for the Nigerian Bitumen deposit. Journal of Petroleum Science and Engineering, 116: 28-35. Lie, K.A. et al., 2012. Open-source MATLAB implementation of consistent discretisations on complex grids. Comput Geosci, 16(2): 297-322. Manchuk, J.G. and Deutsch, C.V., 2013. Optimization of Drainage-Area Configurations To Maximize Recovery From SAGD Operations. Journal of Canadian Petroleum Technology 52(3): 233-242. 23
Rahim, S., Li, Z. and Trivedi, J., 2015. Reservoir Geological Uncertainty Reduction: an Optimization-Based Method Using Multiple Static Measures. Mathematical Geosciences, 47(4): 373-396. Shakhsi-Niaei, M., Iranmanesh, S.H. and Torabi, S.A., 2014. Optimal planning of oil and gas development projects considering long-term production and transmission. Computers & Chemical Engineering, 65(0): 67-80. Tarhan, B., Grossmann, I.E. and Goel, V., 2009. Stochastic programming approach for the planning of offshore oil or gas field infrastructure under decision-dependent uncertainty. Industrial and Engineering Chemistry Research, 48(6): 3078-3097. Tavallali, M.S., Karimi, I.A., Halim, A., Baxendale, D. and Teo, K.M., 2014. Well Placement, Infrastructure Design, Facility Allocation, and Production Planning in Multireservoir Oil Fields with Surface Facility Networks. Industrial & Engineering Chemistry Research, 53(27): 11033-11049. Tavallali, M.S., Karimi, I.A., Teo, K.M., Baxendale, D. and Ayatollahi, S., 2013. Optimal producer well placement and production planning in an oil reservoir. Computers & Chemical Engineering, 55(0): 109-125. Teare, M., Cruickshank, R., Miller, S., Overland, S. and Marsh, R., 2013. Alberta’s Energy Reserves 2013 and Supply/Demand Outlook 2014–2023, Energy Resources Conservation Board, Government of Alberta: Calgary, Alberta. Vanegas, J.W., Cunha, L. and Deutsch, C.V., 2011. Proxy models for fast transfer of static uncertainty to reservoir performance uncertainty. AAPG Memoir(96): 203-215.
24
Table 1. Fixed DA and SP dimensions Parameter
Value (m)
DA length
850
DA width
600
SP length
150
SP width
50
DA–SP spacing
50
Table 2. First stage optimization results Varialbe
Value
θ (rad)
0.7229
∆x (m)
-64.45
∆y (m)
-273.1
Available Bitumen (NCB)
16720
Table 3 First stage optimization results using decomposition scheme 1 Optimization results
Area 1
Area 2
Area 3
θ (rad)
0.6283
0.7362
0.5642
∆x (m)
71.05
25.25
-66.45
∆y (m)
386.2
-169.5
-395.6
Available Bitumen
4056
5646
5970
Total Available Bitumen
15681
25
Table 4. First stage optimization results using decomposition scheme 2 Optimization results
Area 1
Area 2
Area 3
θ (rad)
-0.9675
0.9362
0.6429
∆x (m)
103.55
108.75
-178.45
∆y (m)
-421
226.7
59.4
Available Bitumen
4082
7063
5802
Total Available Bitumen
16947
Table 5. Athabasca reservoir characteristics (Guo et al., 2013) Parameter
Value
Parameter
Value
w (m)
75
a
1.9
k (μm2)
2.5
b
0.66
α (m2 s-1)
7.06×10-7
β
0.264
m
4
γ
0.132
νs (mm2 s-1)
10
ur
0.556
S0
0.8
ud
0.864
ϕ
0.35
vr
2.9
h (m)
26.5
vd
2.125
Table 6. Reservoir characteristics Parameter
Value
DAcell
510
dp (m)
90×10-6
Av (m-1)
6.67×10+4
τ (%)
81
26
Table 7. Design and economic parameters (Giacchetta et al., 2015; JacobsConsultancy, 2012; Lawal, 2014) Parameter
Value
Parameter
Value
DAwell
6
Discount Rate (%)
6
L (m)
850
Tax Rate (%)
30
SOR
2.7
OP (CAN$ m-3)
719.13
3
-1
-3
Basic CPF Capacity (Mm yr )
3.5
FC (M$ m )
3.023×10-7
SAGDlifetime (yr)
25
CCCPF (M$)
2480
Royalty (%)
10
CCDA (M$)
78
Table 8. Solution report for each case study c2
CPU
Optimality
(%)
time (s)
gap (%)
50
3600
1.56
100
804
0.98
150
0.28
0.78
200
0.14
0.03
250
0.02
0.01
300
0.09
0.45
27
Figure 1. SAGD well illustration (left: 3-D view of a SAGD well pair; right; 3-D view of a drainage area with 3 well pairs)
Figure 2. Overview of Optimization Framework 28
Figure 3. Net Continuous Bitumen Calculation
Figure 4. Total NCB Calculation of a DA from aerial view
29
Figure 5. Leased area with surface restriction of the reservoir for well placement
Figure 6. Porosity distribution of the leased area
30
Figure 7. Permeability distribution of the leased area
Figure 8. Net Continuous Bitumen profile
31
Figure 9. Compact DA placement plan
Figure 10. Compact DA placement plan using decomposed area (scheme 1)
32
Figure 11. Compact DA placement plan using decomposed area (scheme 2)
Figure 12. Predicted oil production rates for the Athabasca reservoir
33
Figure 13. Case study 50%–solution
Figure 14. Case study 150%–solution 34
Figure 15. Case study 200%–solution
Figure 16. Optimal NPV at different scenarios for CPF capacity 35
Figure 17. Planning developments on the compact DA placement plan
Figure 18. Oil production rate of DAs over lifetime of the project
36
Highlights
Strategic planning of oil sands development through SAGD process Two-stage optimization framework to plan the development of drainage areas Maximization of the recoverable bitumen through drainage area layout optimization Maximization of the net present value through mixed integer linear optimization
37
PII: DOI: Reference:
S0920-4105(15)30187-X http://dx.doi.org/10.1016/j.petrol.2015.11.023 PETROL3256
To appear in: Journal of Petroleum Science and Engineering Received date: 24 June 2015 Revised date: 17 October 2015 Accepted date: 25 November 2015 Cite this article as: Hossein Shahandeh, Shahed Rahim and Zukui Li, Strategic Optimization of the Oil Sands Development with SAGD: Drainage Area Arrangement and Development Planning, Journal of Petroleum Science and Engineering, http://dx.doi.org/10.1016/j.petrol.2015.11.023 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Strategic Optimization of the Oil Sands Development with SAGD: Drainage Area Arrangement and Development Planning Hossein Shahandeh, Shahed Rahim, Zukui Li
Department of Chemical and Materials Engineering, University of Alberta, Edmonton, AB, T6G2V4 Canada
Abstract The majority of the oil sand deposits in Alberta Canada can only be extracted using in situ methods, mainly the Steam Assisted Gravity Drainage (SAGD). SAGD project consists of Central Processing Facility (CPF) and Surface Pad (SP) situated surface and Drainage Area (DA) consisting of multiple injector and producer well pairs situated subsurface. DA placement is crucial in ensuring that maximum amount of bitumen is extracted. In this work, an optimization framework is developed to plan the development of DAs. The first step aims to obtain a compact DAs arrangement by maximizing the amount of extractable bitumen. In the second step, a Mixed Integer Linear Optimization (MILP) model is developed to arrange the multiperiod development plan of the DAs. The proposed method is applied to a case study with multiple DAs for arrangement and the results demonstrate that the method can effectively generate a good DA layout and an economically optimal development plan that maximizes the Net Present Value (NPV).
Keywords Oil sands, SAGD, reservoir development planning, mixed integer linear optimization
Author to whom correspondence should be addressed. Email: [email protected]; Tel: 1-780-492-1107; Fax: 1-780-492-2881 1
1. Introduction Alberta’s oil sands are the third largest proven global reserves of oil with current estimates of 26.6 billion cubic meters of crude oil (Teare et al., 2013). Most of the oil sands in Alberta are discovered in three geological areas: Athabasca, Peace River and Cold Lake. Among them, Athabasca is the largest in terms of area. About 20% of the bitumen reserves are extractable using surface mining methods. The rest 80% are more than 70 meters below the ground and are extractable using in situ methods. Prominent in situ methods currently used in the oil sands industry are Cyclic Steam Simulation (CSS) and Steam Assisted Gravity Drainage (SAGD), whereas SAGD accounts for the most widely used in situ methods in Alberta. Nowadays, more than 16 SAGD projects are under operation in Alberta and many more are under development. Some examples of the prominent SAGD operations with significant bitumen production are: MacKay River and Firebag operations of Suncor Energy, Foster Creek and Christina Lake operations of Cenovus Energy and Jackfish operations of Devon energy (Teare et al., 2013).
SAGD operations use super-heated steam to decrease the viscosity of the bitumen in the reservoir. As a result, the bitumen can easily flow and be pumped up to the surface. As shown in Figure 1, SAGD wells consist of two horizontal wells which are in parallel. Steam is injected to the upper injector well and the melted bitumen flows to the lower producer well by gravity (Butler, 1991). The Bitumen from the producer well is pumped to the surface and transported to Central Processing Facility (CPF). Except CPF, the major SAGD facilities include Surface Pad (SP) and Drainage Area (DA). SP is the surface facilities from which multiple producer and injector well pairs are drilled. DA includes a set of parallel wells attached to the SP. CPF is a surface facility which processes the oil-water emulsion, recycles the water and produces the super-heated steam required by the injector wells. Steam is supplied from the CPF to each of the SP through pipelines. SAGD well pairs are drilled in areas with high volume of recoverable bitumen, which depends on the porosity and the oil saturation of the reservoir.
Current practices of commercial placement of SAGD wells largely rely on reservoir characterization and engineering judgement. Optimization methods for the placement of SAGDs’ DA and SP are relatively new in academia. Kumar investigated a space packing optimization method for generating a compact and non-overlapping arrangement of DAs (Kumar, 2011). Geometric transformations of global rotation, global translation and column translation were 2
used by the space packing algorithm to maximize the recoverable bitumen. Restrictions such as non-placement of SP over surface restrictions or non-placement of DA in thief zones were considered by incorporating a penalty function in the objective function. Manchuk and Deutsch used an adaptive grid search algorithm to determine the optimal arrangement of DAs (Manchuk and Deutsch, 2013). The optimization algorithm determined the positions and orientations of the SPs and DAs over a reservoir area to economically maximize the recoverable bitumen.
During the lifecycle of an oilfield project, development is the most critical phase due to intensive investments required for capital and operating costs. Determining whether it is economic to develop an oilfield or not and finding the best schedule for field development are the main concerns in this field of study. According to the timeframe, optimal decision making can be categorized as operational, tactical and strategic ones. On the other hand, exploration, development, production and transmission represent the major decision scopes (Shakhsi-Niaei et al., 2014).
Using simplifying assumptions such as linear performance of reservoirs versus time, optimization of the oilfield investment and operation planning was first tackled with Linear Programming (LP) (Aronofsky and Williams, 1962). Simultaneous optimization of locationallocation of wells, facility operation scheduling, and production rates during different time periods was then modeled as a Mixed Integer Linear Optimization (MILP) problem, but still based on linear performance of reservoirs (Frair, 1973). There have been many contributions in this field of study then after, and they can be mainly categorized into offshore and onshore oilfields’ development. Emphasizing on providing a detailed and accurate geological model of reservoirs, Behrenbruch outlined a systematic procedure for feasibility study of developing offshore petroleum fields (Behrenbruch, 1993). Tarhan et al. developed a multistage stochastic programming model for the planning of offshore oil/gas field infrastructure under uncertainty (Tarhan et al., 2009). Recently, Gupta and Grossmann proposed Mixed Integer Nonlinear Programming (MINLP) formulations for offshore oilfield development planning (Gupta and Grossmann, 2012; Gupta and Grossmann, 2014). In the former work, facility installation, well drilling and production decisions considering oil, water and gas flows profiles were simultaneously optimized when the problem was reduced to a MILP (Gupta and Grossmann, 2012). In the later work, a general multistage stochastic programming was proposed including
3
endogenous uncertainties (either fiscal rules or uncertainty in the field parameters) (Gupta and Grossmann, 2014).
Ierapetritou et al. proposed a novel decomposition-based approach for field development problem based on field decomposition and a quality cut-off criterion based iterative scheme (Ierapetritou et al., 1999). The optimization problem could be solved to optimality for each of the decomposed subproblem when it was formulated as a MILP problem. Kosmidis et al. presented a novel MINLP model considering nonlinear reservoir behaviour, multiphase flow in wells and constraints from the surface facilities (Kosmidis et al., 2005), by inclusion of piecewise linear approximations of each well model and using an outer approximation based algorithm. Tavallali et al. tried to apply a more realistic and detailed reservoir model and consequently efficient solution algorithm in two consecutive studies (Tavallali et al., 2014; Tavallali et al., 2013). The major contribution of the first work was the consideration of rigorous subsurface flow dynamics (Tavallali et al., 2013). They also introduced a unified model integrating subsurface, wells, and surface levels of an upstream production project. In the second work, a model was developed for location-allocation problem of well-drilling and infrastructure installation when the oilfield has more than one reservoirs with a shared surface processing network (Tavallali et al., 2014). Both the interactions within the surface network facilities and the dynamic states of the subsurface environment were incorporated into the model.
As the conventional oil resources become unavailable or uneconomic for extraction, bitumen has become more and more favorable. As an in situ recovery method, SAGD can be an excellent case study in the scope of strategic development optimization since its facility and recovery mechanism are very different from the conventional methods. Continuously injecting steam, having rising and depletion stages during bitumen recovery, and distinctive geological characteristics of each DA make this process a unique one. With this motivation, an optimization framework for determining the optimal placement and development order of SAGDs’ DA and SP facilities is proposed in this study. The objective of the DA arrangement optimization is to determine DA and SP locations which would provide maximum accessibility of the available bitumen resources. The objective of the development plan optimization is to maximize the Net Present Value (NPV). The optimal SAGDs’ DA and SP arrangement and development plan are obtained through a two-stage optimization framework. In the proposed method, the derivative free optimization solver NOMAD is used in Matlab to perform DA configuration optimization, 4
the mixed integer linear optimization solver CPLEX is used in GAMS to perform the development planning optimization.
The rest of the paper is organized as follows. The problem statement and the workflow for the proposed framework are explained in section 2. The proposed DA arrangement optimization method is presented in Section 3 followed with a case study. Section 4 presents the optimization model and results for reservoir development planning. The paper is concluded in section 5.
2. Problem statement
There are two objectives related to the strategic planning of the SAGD development. The first objective is to obtain a compact arrangement of DAs within the reservoir lease area to achieve maximum bitumen recovery. The second objective is to generate a plan for developing the DAs (i.e., determine the starting and ending time within the whole lifetime of the project) in order to reach the optimal NPV. A two-stage optimization method is proposed to achieve the two objectives. An overview of the two-stage optimization method is given in Figure 2 and it is summarized as follows.
The input to the first stage optimization includes the reservoir grid, geological data of the reservoir, any surface restrictions on the reservoir grid and the fixed dimensions of each DA and SP in the compact arrangement. A Matlab program is developed to calculate the total available bitumen with a given DA arrangement. A derivative free optimization method is used for the first stage optimization. The decision variables include the orientation and position of the arrangement of the DAs. The orientation is represented by the angle of the compact DA arrangement from the horizontal x direction. The position is represented by the translation of the compact DA arrangement along both the x and y direction. In the optimization model, a compact and non-overlapping arrangement of DA is considered so as to maximize the accessibility of the available resources. The DAs and SPs are in rectangle shape and two possible locations of SP are considered for each DA. The DA and SP dimensions are considered to be fixed in this study and as a result, fixed numbers of parallel injector/producer wells per DA are considered. The complexity of the optimization model is significantly reduced by fixing those parameters. The optimizer runs till the stop criteria of the optimizer (e.g., maximum number of iterations) is 5
satisfied. The results of the optimizer are the optimal orientation and location of the compact set of DAs. One or two candidate SP locations will be available for each DA in the compact arrangement depending on the surface restrictions and lease boundary.
The output of the first stage optimization provides input to the second stage optimization for generating the DA development plan. Oil production flowrates of all DAs are calculated in Matlab based on rigorous predictive model during DAs' lifetime, and they are used in the second stage optimization. A MILP optimization model is proposed for the second stage of the DA development planning. The MILP optimizer aims to identify the optimal DA development order that maximizes the NPV. This enables an optimal selection of DA for development in different time period. Different scenarios are also considered for the capacity of CPF in the second stage of optimization in order to identify the best capacity and capital investment for the CPF.
3. Drainage area arrangement optimization
Bitumen deposits are found in porous mediums underground between impermeable rocks. The top surface of the bitumen deposit below the impermeable section of rocks is defined as the Top Continuous Bitumen (TCB) and the bottom surface of the bitumen deposit above the impermeable section of rock is defined as the Bottom Continuous Bitumen (BCB). The BCB and TCB of a bitumen deposit are determined from data obtained from well logs or exploratory drills. Net Continuous Bitumen (NCB) is the total thickness between the TCB and BCB that meets certain minimum criteria for reservoir quality (Kumar, 2011) and it is calculated as the number of net cells which are above a threshold porosity and permeability value. For each DA, the available bitumen is the quantity of bitumen available for production along all the wells within the DA. In this study, the NCB value is used to quantify the available bitumen. Figure 3 illustrates how NCB is calculated along the elevation of the reservoir, where the black portions of the reservoir have lower porosity and permeability value and therefore they are not included in the determination of NCB. For the ith DA, the available bitumen is calculated by the following equation
Riavailable
NCB
cDAi
c
arc So
(1)
6
where Riavailable is the available bitumen, So is the oil saturation, NCBc is the net continuous bitumen for cell c. arc is the fraction of area of cell c inside the ith DA and its calculation is illustrated in Figure 4. Constant oil saturation So is assumed for the whole reservoir in this work. As a result, the total NCB value calculated as
cDAi
NCBc arc is a good indication of the
amount of bitumen available for extraction by the SAGD wells. Since the oil saturation So appears as a constant coefficient in (1), its value does not affect the DA configuration. To optimize the DA configuration and maximize the total available bitumen, a Matlab program was developed which calculated the total available bitumen as a function of the orientation and position of the compact set of DAs. The input of the program includes:
Grid size with the lease area boundaries and surface restriction on the grid,
Geological data for each cell in the grid,
DA dimensions, SP dimensions, DA-SP spacing,
Orientation (θ) and positions (∆x, ∆y) of the compact arrangement of DAs.
For each DA arrangement plan defined by (θ, ∆x, ∆y), the total available bitumen is evaluated based on the following principles:
A DA is considered valid only if its area is 80% or more inside the lease boundary. For a rectangle shape DA, we assume it is worth to be developed only when the majority of its area is within the lease boundary. Otherwise, the major part of the DA will be outside the lease boundary such that its development with horizontal well is not practical.
A DA is considered valid only if its total available bitumen is above a certain threshold value. This assumption is made to avoid placing a DA in the region that has very low available bitumen and is not economical to be developed. This threshold is determined using the following method. Initially, a DA configuration with zero angle and zero translation is evaluated to find the distribution of the available bitumen for all the DAs. Based on the distribution, the threshold value for total available bitumen is set as the number such that 10% of the values lie below it.
A DA is considered valid only if it has at least one possible location for placing a SP. For a feasible SP location, the SP had to be placed inside the lease boundary and could not be placed on top of any surface restrictions (e.g., rivers, lakes, roads or conserved forest areas).
The output of the Matlab program is the total available bitumen from the given DA arrangement. 7
The first stage in the optimization framework generates a compact arrangement of DA. In this step, the objective function for the compact DA arrangement optimization problem is designed as maximizing the total available bitumen, which can be represented as
max Riavailable ( , x, y )
(2)
i
where Riavailable is the available bitumen for the ith DA. However, there is no explicit equation for the above objective function since it is evaluated through a Matlab program. Derivative free optimization technique is used in this work which calls on the developed Matlab program for total available bitumen calculation. Specifically, the derivative free optimization solver NOMAD is used, which implements the Mesh Adaptive Direct Search (MADS) algorithm for constrained blackbox functions. The MADS algorithm is an extension of pattern search method for nonlinear constrained optimization problems (Audet et al., 2009). During this optimization, upper and lower bounds were applied for the decision variables: orientation, translation along x direction and translation along y direction (by assuming the left bottom corner of the reservoir grid is initially located at the origin with x direction parallel to the horizontal axis), as given by Equations 3–5 respectively.
(3)
x x x
(4)
y y y
(5)
where is radian, x and y are twice the length and width of the reservoir, respectively.
Case Study
A case study with a three dimensional reservoir grid is used to demonstrate the practicality of the optimization method used for SAGD DA placement. A relatively large reservoir model with a grid size of 55 × 62 × 10 cells with cell dimensions of 100m × 100m × 5m is considered. The geological data considered for the case study are the porosity and permeability of the reservoir. Porosity and permeability were generated for every cell of the grid using a built-in function GaussianField in Matlab Reservoir Simulation Toolbox (Lie et al., 2012). The porosity values were in the range of 0.1 to 0.5. Permeability was calculated from the porosity data using Carmen-Kozeny relationship. Furthermore, in this case study an irregular shaped leased 8
boundary was used to reflect the common situation of oil sands projects. In addition, surface restrictions were created within the lease boundary to reflect surface restrictions for placing SP. The aerial view of the lease boundary and the surface restriction is given in Figure 5. In Figure 5, the leased area is represented as the yellow section, surface restriction is represented by the blue area and the white parts represent the area outside the lease boundary. Notice that the surface restriction area is virtual in this case study, and it has been divided into convex polygons for use in the Matlab code of NCB evaluation. Every SP is forced to have no interaction with any of those polygons. However, DA can be placed under those surface restrictions since it is underground.
The porosity and permeability distribution of the three dimensional grid is given in Figure 6 and 7 respectively. The dark blue regions in those figures represent the unleased area. Geological data for the unleased area is considered to be zero to simplify the recoverable bitumen calculation.
The geological data is used to calculate the NCB of the reservoir. In this case study, any cell with a porosity and permeability above a threshold value of 0.25 and 1.1×10-13 m2 respectively, were considered to be a net cell. A threshold value for the total bitumen is applied for each DA to be counted in the objective function (2). It is determined as 300 in this case study using the aforementioned method. The aerial view of the NCB of the reservoir grid is given in Figure 8.
In order to provide SAGD well arrangement plan, a list of parameters specific to all the DAs and SPs in the compact arrangement were fixed. The fixed dimensions of the DA and SP are given in Table 1. Based on the setting, each DA has six pairs of injector and producer wells. Two possible locations of SP with respect to a specific DA were considered with a fixed DA-SP spacing.
Overall placement plan The proposed optimization method was applied to the reservoir and it provides a compact arrangement of DAs. The results of orientation and translation along the x and y direction for the compact DA layout is given in Table 2. The total available bitumen in Table 2 is the sum of the available bitumen values of all the DAs in the compact arrangement.
9
The compact arrangement of all the DAs from the first stage optimization is given in Figure 9. The numbers on each of the DA quantifies the amount of available bitumen (NCB value) from that particular DA. The following observations can be made from Figure 9:
No DA is placed outside of the leased area
No DA is placed on areas which would result in a significant fraction of DA area to be outside the lease boundary
No SP is placed outside of the leased area
No SP is placed on any of the surface restrictions
It is important to note that due to surface restriction and lease boundaries, some DAs have only one possible location for the placement of SP.
Decomposition based placement The first stage optimization method was further studied by decomposing and breaking down the lease area into subareas. The DA configuration optimization was applied to each of those subareas to determine the optimal orientation and position for the compact set of DAs. Two lease area decomposition schemes are investigated here. The first scheme is shown in Figure 10, where the red dotted lines divide the lease area into 3 subareas. Table 3 lists the optimal orientation and position of the different DA arrangement for 3 subareas. The total available bitumen of each of the subareas and the sum of the total available bitumen from all the subareas are also reported in Table 3. The results show that no DA or SP is placed outside the lease boundary and no SP is placed on any of the surface restricted areas. Furthermore, a DA was not placed if the DA did not have any available position for placement of at least one SP. For the above decomposition scheme, the total available bitumen from the 3 divided areas (15,681) is lower than the total available bitumen from the undivided compact DA arrangement plan of Figure 9 (16,720).
The second decomposition scheme is shown in Figure 11. The DA arrangement optimization was then performed on the 3 decomposed sections and the total available bitumen from the 3 different sections are reported in Table 4. In this case, the total available bitumen from summation of the total available bitumen from all the different areas (16,947) is higher than the total available bitumen from the DA arrangement obtained from the undivided lease area (16,720). The above study shows that decomposing a large lease area into small subareas and then applying the compact DA arrangement optimization method separately to those subareas may result in higher total available bitumen. However, it requires testing different decomposition 10
scheme and finding the best way to divide the lease area. This can be deemed as another optimization problem for seeking the global optimal DA configuration. Consider the complexity of existing optimization problem, this study is not included in this work.
4. Reservoir development planning
The objective of the second stage optimization of the proposed framework is to determine the DA development plan based on the compact arrangement solution obtained from the first stage. In this section, the oil production rate profile is first calculated through simulation of nonlinear production models for each DA according to its geological specifications. Parameters used in the calculation are provided and explained. A MILP problem is then formulated for planning the DA development with the objective of maximizing the NPV. The DA layout shown in Figure 9 was used as the basis for the case study of this section.
4.1 Well production model
The crucial part of the second stage optimization is a model for oil production rate through SAGD process based on its geological characteristics. Some models have been developed for this purpose in which the SAGD production is divided into rising steam chamber stage and lateralspreading steam chamber stage (Akin, 2006; Guo et al., 2013; Vanegas et al., 2011), and the oil production rates in the two different stages are calculated as
akgSoh t 2 kg / m s qd 2 2 2 m s w bSoh 1/2
3/2
4 qr 3
9 42
2/3
(6)
2/3
So
1/3
kg 1/3 t m s
(7)
where qd and qr (m2 s-1) are oil production rate per meter of well length in depletion and rising stage respectively, k (m2) is the effective reservoir permeability, g (m s-2) is the acceleration due to gravity, α (m2 s-1) is the reservoir thermal diffusivity, ϕ is the porosity, ΔSo is the displaceable oil saturation, h (m) is the effective drainage height, m is a viscosity-temperature correlation parameter, νs (m2 s-1) is the oil kinematic viscosity at the steam temperature, t (s) is the time, w (m) is the half well spacing, γ is the shape factor of the rising steam chamber, and a, b, β are
11
empirical constants. Among the mentioned characteristics, ΔSo is the only one varying with time (Guo et al., 2013);
h So So 0.43 s kgt
0.4
(8)
where So is initial oil saturation of the reservoir. To predict the oil production rate, the depletion model and rising model are combined together with weighting factors (Guo et al., 2013):
qtoil yd qd yr qr
(9)
yd
1 1 10ud ( vd t )
(10)
yr
1 1 10ur ( vr t )
(11)
where qtoil is the total oil production rate, yd is the weighting factor of depletion model increasing from 0 to 1, yr is the weighting factor of rising model decreasing from 1 to 0, t (year) is the time, and ud, vd, ur, vr are constants. In the above equations, there are many nonlinear terms (especially for depletion stage) that will cause difficulty for optimization. Therefore, they were not directly used in the optimization model for oil production rate calculation. Instead, oil production rates were first calculated versus time for each DA, and the rates were then utilized as parameters to describe the oil production rate at each year. In this way, the optimization problem could be simplified and then solved with much less computational efforts.
Athabasca reservoir characteristics adapted from reference (Guo et al., 2013) were used to validate the presented model for oil production rate and also demonstrate trend of qd, qr, and qtoil versus time. Related parameters are provided in Table 5. Figure 12 shows the rising stage is the dominant stage at the initial stage and the depletion stage is the dominant one at the end, while they are both partially effective in the middle.
For the case study of this work, oil production rates of each individual DA can be predicted according to its reservoir characteristics. Average porosity, permeability, and initial oil saturation
12
are the parameters which can be defined for each DA based on the following equations (Rahim et al., 2015);
0.25 NCB DAcell
Av k
(12)
6 dp
(13)
1 3 2 Av 2 1 2
100 2.25 k Sw
(14)
2
(15)
So 1 Sw
(16)
where DAcell is the number of cells in each DA, Av (m-1) is the specific area of spherical shape grains, dp (m) is the average diameter of grains, τ is tortuosity, and Sw is water saturation. Constants are reported in Table 6. It should be highlighted here that each DA was assumed to have six pairs of injector and producer wells.
4.2 Planning optimization model
Based on the above well production model, the DA development planning problem can be defined next. In the model, indices i, t and k correspond to DA, time period within the lifetime of the whole project, and time period within the lifetime of a DA, respectively. The total number of DAs was 43 and the total lifetime of the project was set as 25 years. The total lifetime of each DA, DAiLife , was obtained from Matlab calculation and it varies from 6 to 10 years depending on its geological property. Note that 10 years is an upper bound of the life time for all the DAs based on the production rate calculation.
In the optimization model, the starting point ( tiStart ) and ending point ( tiEnd ) of each DA are the main decision variables. It is assumed that decisions can only be made at the beginning of each year and any operation of a DA can only be changed at that time. Starting point is defined by following equation: tiStart xi ,t t i
(17)
t
13
where xi,t is a binary variable which takes value 1 if ith DA starts its operation at time t, and value 0 otherwise.
To make sure that there is only one starting point for each DA, the next constraint is applied. If the summation of binary variables is equal to zero, it means that DA should not be operated at all during the lifetime of the project.
x
i ,t
1 i
(18)
t
The lifetime of each DA, tiLife , is a parameter varying from 6 to 10 based on its geological property. The ending point of a DA can be then determined as
tiEnd tiStart tiLife
x
i ,t
i
(19)
t
If being under operation of a DA is uneconomical, corresponding starting and ending times will be both equal to zero. There should be some mathematical logics between the oil production rate qioil,k of ith DA at kth year of its lifetime and the oil production rate of the same well in 25 years basis of the whole project f i oil ,t . In the other words, the lifetime of each DA can be shifted through the project horizon. This can be achieved by introducing a new binary variable yi ,k ,t . For each DA i, therefore, there exists a matrix of binary variables with k rows and t columns which can only have a value of one in a diagonal form. To determine the value of one in the first row of each matrix a big-M formulation can be applied and the oil production rate qioil,k would be equal to f i oil ,t if only t is equal to tiStart k 1 . The following equation was applied to reflect this logic: 1 yi ,k ,t M t k 1 tiStart 1 yi ,k ,t M t k 1 i, t, k 1
(20)
Equation 20 is a “big-M” constraint. It is active when yi,k,t is equal to 1, otherwise, it is redundant. Based on the inequality, M should be bigger than the largest possible value of
tiStart (t k 1) and tiStart (t k 1) , which would be 27 (=19 – 1 + 10 – 1) according to the variables’ domain: a well starting lifetime cannot be later than year 19th and its longest lifetime is 10 years. While any number that is bigger than 27 should be a valid option for M, we choose 28 as the value to keep the constraint tight.
14
For the rest of the rows, the remaining y binary variables can only be equal to one in a diagonal line,
yi ,k ,t yi ,k 1,t 1 i, t 1, k 1, k c1
(21)
As mentioned earlier, DAs have different lifetime varying from 6 to 10 years, but all the y matrices have 10 rows for the sake of mathematical consistency. As a results, 10th, 9th–10th, 8th– 10th, and 7th–10th rows of y matrices must be fixed to zero when the lifetime of corresponding DA is respectively 9, 8, 7, 6 years to avoid any incorrect mathematical logic. The c1 constant number in equation 21 is determined based on this concept. Moreover, there are upper and lower triangles in y matrices which also can set to zero due to diagonal form of binary variable equal to one. In this way, the number of binary variables can be reduced by more than 40%.
With the binary variables yi ,k ,t , the following equation can be used to shift any DA with different lifetime through the project horizon for calculating the oil production rate
yi,k ,t qioil,k k
fi oil ,t L DAwell
i, t
(22)
where L is the length of each well and DAwell is the number of producer wells under operation for each DA. Even though the left hand side of the last equation has a sigma sign, only one of the multiplication will be nonzero due to the way that yi,k,t has been defined. Note that qioil,k , which is calculated through equations 6–11, has the unit of (m3 m-1) and it only corresponds to the 3 flowrate of single pair of wells per its length. On the other hand, unit of f i oil ,t is (m ) and it shows
the whole flowrate of each DA including its six pairs of wells. The above three equations work properly if only the following one is also incorporated. Otherwise, all the yi,k,t would be equal to zero.
y
i ,k ,t
k
tiLife
i
(23)
t
The rest of the model includes economic model of the process which is based on NPV. It was assumed that the main utility consumption of oil production is related to steam. Required steam flowrate and its cost can be estimated as; fi ,steam fi oil t ,t SOR i , t
(24)
TFCt FC fi ,steam t t
(25)
i
15
where f i ,steam is required flowrate of steam for ith DA at tth time period, SOR is a constant number t used for Steam to Oil Ratio, TFC is Total Fuel Cost (M$ yr-1) of the tth year , and FC is Fuel Cost (M$ m-3).
Capital cost of process is composed of CCCPF (M$), initial investment for CPF which is mandatory, and CCDA (M$), related to investment for each DA installation. The value for CCCPF was estimated from the reference (Giacchetta et al., 2015) and a CPF capacity value of 3.5×106 m3 yr-1 was used as the base value, and all the 50%, 100%, 150%, 200%, 250%, and 300% of it were studied as possible scenarios. A proportional relation was implemented between the capital cost of CPF and its capacity. It was also assumed that each DA has 6 pairs of injector and producer wells and capital cost of each of them was 78 M$ (JacobsConsultancy, 2012).
TCCt CCCPF CCDA xi ,t
t 1
(26)
i
TCCt CCDA xi ,t
t 1
(27)
i
Equation 26 is only valid for the first year of project, and equation 27 is valid for the remaining ones. The only difference between equations 26 and 27 is that capital cost of CPF should be just taken into account at the first year. In the following equations, TOMC is Total Operating and Maintenance Cost (M$ yr-1), TAC is Total Annual Cost (M$ yr-1), TAR is Total Annual Revenue (M$ yr-1), OI is Operative Income (M$ yr-1), they can be calculated as (Giacchetta et al., 2015; Lawal, 2014);
TOMC
Discount rate TCCt t
(28)
SAGDlifetime
TACt TFCt TOMC t
(29)
TARt 1 Royalty OP qi ,t t
(30)
i
OI t TARt - TACt t
(31)
where OP is oil price (CAN$ m-3). If a company wants to extract oil sands in an area which is predefined as oil sands area by the government of Alberta, they have to pay a percent of gross revenues as royalty. This number has scaled value according to the price of West Texas Intermediate.
16
Tax can also be calculated using following equation (Giacchetta et al., 2015; Lawal, 2014);
Taxt TOI t Tax rate t
(32)
TOI t OI t Akt
(33)
Akt
t
xt CCCPF CCDA i ,Life SAGDlifetime i DAi
Akt CCDA i
xi ,t DAiLife
t 1
t 1
(34)
(35)
where TOI is Total Operating Income and Ak is amortization of investment.
NPV is calculated based on NI (Net Income) and DCF (Discount Cash Flow). They can be calculated as (Giacchetta et al., 2015; Lawal, 2014):
NI t OI t Taxt t
DCFt
(36)
NI t
1 Discount rate
t
t
(37)
NPV DCFt TCCt t SAGDlifetime t
(38)
t
Two additional constraints were imposed in this optimization problem. The first one was on the total amount of DAs installation’s investment at each year which should be equal or less than 300 M$.
TACt TCCt 300 t
(39)
Furthermore, total oil production at each year must be equal or less than the capacity of CPF. When a capacity is selected for a plant during design phase, processing of raw materials more than that fixed number would be impossible unless using the existing equipment more effectively or purchasing new equipment. While increasing the capacity has not been considered here as an option, the following constraint could be used as the upper limit of CPF capacity:
f
oil i ,t
3.5 106 c2 t
(40)
i
where c2 constant is a coefficient to consider different scenarios of 50%, 100%, 150%, 200%, 250% and 300% of the nominal capacity (3.5×106 m3/year) of CPF.
Finally, the optimization model can be expressed as 17
max NPV s.t. eqs(17 40)
(41)
The optimization problem (41) is a MILP problem. With all the required constants provided in Table 7, the model was solved in GAMS using the CPLEX solver. The optimization stopping criteria were set as one hour computational time and relative optimality gap of 1%. The solution for the base capacity of CPF (i.e., c2=1) was found in less than 804 seconds with a desktop computer (single core of Intel® i5-4590 @ 3.30 GHz, 8 GB RAM). The optimal NPV value at the end of the project lifetime was 12559.6 M$.
4.3 Result and discussion
As mentioned in previous section, 6 different scenarios for the capacity of CPF were taken into account to investigate the design characteristic of the project. Computational statistics for each case study are available at Table 8. To have a detailed study, the results of three capacity scenarios (50%, 150% and 200%) are illustrated through Figures 13–15. Each of them consists of a development plan (top left), a profile of the number of operating DAs (top right), NPV (bottom left) and total oil production (bottom right) over 25 years of project lifetime. Note that the segments in the development plan figure represent the operating years of each DA.
In this study, all the 43 DAs shown in Figure 9 were first sorted based on their NCB (bitumen) values in an ascending order. Therefore, the DA number 1 has the lowest NCB and the DA number 43 has the highest NCB. As the Figures 13-15 show, the DAs with higher NCBs are more favorable for being developed and under operation at the early stage of the project lifetime. When the capacity of CPF is small, many DAs with small NCBs are not even applied at all. On the other hand, when the CPF has a large capacity, all the DAs tend to be developed while the priority of operation is still on the DAs with high NCB values. If the DA’s investment constraint and CPF capacity constraint had not been imposed, all the DAs would have been developed and operated since the first year. The reason is that keeping the operation running for a longer time means higher amortization and operating cost. Activation of the investment limit constraint and CPF capacity constraint forces the project to take longer time of operation, and only three DAs could be installed during each year. Another interesting fact about these results is that the DA’s investment constraint was active at c2=50%, 150%, 200%, 18
while the constraint associated with the capacity of CPF was completely active at c2=50%, partially active at c2=150%, and completely inactive at c2=200%. Due to CPF investment in the first year, NPV starts to grow from a negative value. Afterwards, it becomes higher and higher versus operating year since the capital cost of DAs’ installation is not that much compared to the oil production revenue. In the last years, NPV does not change significantly since the production rate decreases more and more and there is not any investment for a new DA in the last five years.
The larger the CPF capacity, the higher number of DAs could be operated and also the higher total oil production could be achieved. Nevertheless, larger CPF requires higher capital cost investment which will lead to a lower NPV. As demonstrated in Figure 16, there exists an optimal point for the capacity of CPF which happens at 5.25 Mm3 yr-1 (150% of the base capacity).
Since the optimal scenario was found at this point, it would be beneficial to provide more information about it such as planning development of the case study on compact DA placement plan (Figure 17) and flowrate of each individual DA during project lifetime (Figure 18). Figure 17, which is based on the DA layout in Figure 9, shows different DAs and corresponding NCB (black number) and lifetime (yellow number) to them. For example, the DA at the top left corner has the NCB of 309 and is under operation from year 14th to 24th. It worth to mention that darker rectangular represent DAs with higher NCB and lighter ones represent DAs with lower NCB.
Oil production of each DA versus operating year is also provided in Figures 18. As it can be observed, DAs with better performance are applied for operation at beginning years. DAs’ oil production profiles are represented in Figure 18 and darker lines belong to DAs with larger ID (i.e., higher NCB values).
5. Conclusion
A novel optimization framework for strategic planning of in situ oil sands production has been proposed in this paper. The proposed method determines the compact arrangement of SAGD’s DAs and the optimal development order for those DAs. The first stage of the proposed 19
optimization method determines the optimal orientation and position of a compact set of all the DAs and SPs that can be packed within the lease area of the reservoir so that the total extractable bitumen is maximized. Surface restrictions and irregular shaped lease areas are taken into consideration in the placement of DAs and SPs. The SAGD arrangement optimization model produced a compact DA arrangement plan with high recoverable bitumen resource. The second stage of the proposed optimization method determines the optimal DA development plan for the lifetime of the project. This model generates a plan that maximizes the NPV under the investment and CPF capacity restrictions. Different scenarios were also studied for capacity of CPF to show the importance of this designing factor. The optimization method was applied to a case study with realistic reservoir size to determine a DA placement and development plan that maximizes the NPV from the SAGD project.
There are several ways to further improve the proposed method. First, the proposed DA configuration optimization problem is based on fixed dimensions of each DA. To achieve global optimal placement, non-uniform size of each DA should be modeled although it leads to a more complex optimization problem since the number of decision variables significantly increase. In addition, developing a optimization model for identifying the best possible decomposition of lease area is also a way to further improve the solution. Finally, in future work, the SAGDs’ DA arrangement optimization method will be studied by incorporating geological uncertainty (i.e., multiple geological realizations). Long term uncertainty in crude oil demand and price will also be considered to study the sensitivity and robustness of the reservoir development plan.
Acknowledgements The authors would like to acknowledge the financial support from Natural Sciences and Engineering Resource Council (NSERC) of Canada Discovery Grant Program.
20
Nomenclature Indices i
drainage area (DA)
c
cell
k
time period within the lifetime of a DA
t
time period within the lifetime of the whole project
Parameters a, b, β
empirical constants in oil production rate equation in depletion stage
Av
specific area of spherical shape grains (m-1)
CCCPF
initial investment for CPF (M$)
CCDA
investment for each DA installation (M$)
DAcell
number of cell in each DA
DAlife
total lifetime of each DA
DAwell
number of producer wells under operation for each DA
dp
average diameter of grains (m)
FC
Fuel Cost (M$ m-3)
g
acceleration due to gravity (m s-2)
h
effective drainage height (m)
k
effective reservoir permeability (m2)
L
length of each well (m)
m
viscosity-temperature correlation parameter
M
big-M constant
OP
Oil price (CAN$ m-3)
So
initial oil saturation of the reservoir
SOR
Steam to Oil Ratio
Sw
water saturation
ud, vd
empirical constants for weighting factor equation of depletion model
ur, vr
empirical constants for weighting factor equation of rising model
w
half well spacing (m)
α
reservoir thermal diffusivity (m2 s-1)
γ
shape factor of the rising steam chamber
νs
oil kinematic viscosity at the steam temperature (m2 s-1)
τ
tortuosity
porosity
Continuous variables Ak
amortization of investment 21
ar
fraction of area of cell c inside the ith DA
DCF
Discount Cash Flow
foil
oil production rate in 25 years basis of the whole project (m3 s-1)
fsteam
required flowrate of steam
NCB
Net Continuous Bitumen
NI
Net Income
OI
Operative Income (M$ yr-1)
qd
oil production rate per meter of well length in depletion stage (m2 s-1)
qoil
oil production rate in lifetime years basis of each DA (m2 s-1)
qr
oil production rate per meter of well length in rising stage (m2 s-1)
qtoil
total oil production rate (m2 s-1)
Ravailable
available bitumen for each DA
t
time (s)
TAC
Total Annual Cost (M$ yr-1)
TAR
Total Annual Revenue (M$ yr-1)
tend
ending point
TFC
Total Fuel Cost (M$ yr-1)
tlife
lifetime of each DA
TOI
Total Operating Income
TOMC
Total Operating and Maintenance Cost (M$ yr-1)
tstart
starting point
yd
weighting factor of depletion model
yr
weighting factor of rising model
ΔSo
displaceable oil saturation
∆x
horizontal position of the compact arrangement of all DAs
∆y
vertical position of the compact arrangement of all DAs
θ
orientation of the compact arrangement of all DAs
Binary variables xi,t
whether the operation of ith DA starts at time t or not
yi,k,t
whether the operation of ith DA starts at time k (DA lifetime) and time t (project lifetime) or not
22
References
Akin, S., 2006. Mathematical modeling of steam-assisted gravity drainage. Computers and Geosciences, 32(2): 240-246. Aronofsky, J.S. and Williams, A.C., 1962. The Use of Linear Programming and Mathematical Models in Underground Oil Production. Management Science, 8(4): 394-407. Audet, C., Digabel, S.L. and Tribes, C., 2009. NOMAD user guide, Les Cahiers du GERAD. Behrenbruch, P., 1993. Offshore oilfield development planning. JPT, Journal of Petroleum Technology, 45(8): 735-743. Butler, R.M., 1991. Thermal recovery of oil and bitumen. Prentice Hall Inc., Englewood Cliffs, NJ. Frair, L.C., 1973. Economic Optimization of Offshore Oilfield Development, University of Oklahoma, Tulsa, OK. Giacchetta, G., Leporini, M. and Marchetti, B., 2015. Economic and environmental analysis of a Steam Assisted Gravity Drainage (SAGD) facility for oil recovery from Canadian oil sands. Applied Energy, 142: 1-9. Guo, J., Zan, C., Ma, D.S. and Shi, L., 2013. Oil production rate predictions for steam assisted gravity drainage based on high-pressure experiments. Science China Technological Sciences, 56(2): 324-334. Gupta, V. and Grossmann, I.E., 2012. An efficient multiperiod MINLP model for optimal planning of offshore oil and gas field infrastructure. Industrial and Engineering Chemistry Research, 51(19): 6823-6840. Gupta, V. and Grossmann, I.E., 2014. Multistage stochastic programming approach for offshore oilfield infrastructure planning under production sharing agreements and endogenous uncertainties. Journal of Petroleum Science and Engineering, 124: 180-197. Ierapetritou, M.G., Floudas, C.A., Vasantharajan, S. and Cullick, A.S., 1999. Optimal location of vertical wells: Decomposition approach. AIChE Journal, 45(4): 844-859. JacobsConsultancy, 2012. Assessment of Innovative Applications of Electricity for Oil Sands Development Phase 1 Report. www.ptac.org/attachments/1425/download. Kosmidis, V.D., Perkins, J.D. and Pistikopoulos, E.N., 2005. A mixed integer optimization formulation for the well scheduling problem on petroleum fields. Computers & Chemical Engineering, 29(7): 1523-1541. Kumar, A., 2011. Optimal drainage area and surface pad positioning for SAGD development, University of Alberta, Canada. Lawal, K.A., 2014. Economics of steam-assisted gravity drainage for the Nigerian Bitumen deposit. Journal of Petroleum Science and Engineering, 116: 28-35. Lie, K.A. et al., 2012. Open-source MATLAB implementation of consistent discretisations on complex grids. Comput Geosci, 16(2): 297-322. Manchuk, J.G. and Deutsch, C.V., 2013. Optimization of Drainage-Area Configurations To Maximize Recovery From SAGD Operations. Journal of Canadian Petroleum Technology 52(3): 233-242. 23
Rahim, S., Li, Z. and Trivedi, J., 2015. Reservoir Geological Uncertainty Reduction: an Optimization-Based Method Using Multiple Static Measures. Mathematical Geosciences, 47(4): 373-396. Shakhsi-Niaei, M., Iranmanesh, S.H. and Torabi, S.A., 2014. Optimal planning of oil and gas development projects considering long-term production and transmission. Computers & Chemical Engineering, 65(0): 67-80. Tarhan, B., Grossmann, I.E. and Goel, V., 2009. Stochastic programming approach for the planning of offshore oil or gas field infrastructure under decision-dependent uncertainty. Industrial and Engineering Chemistry Research, 48(6): 3078-3097. Tavallali, M.S., Karimi, I.A., Halim, A., Baxendale, D. and Teo, K.M., 2014. Well Placement, Infrastructure Design, Facility Allocation, and Production Planning in Multireservoir Oil Fields with Surface Facility Networks. Industrial & Engineering Chemistry Research, 53(27): 11033-11049. Tavallali, M.S., Karimi, I.A., Teo, K.M., Baxendale, D. and Ayatollahi, S., 2013. Optimal producer well placement and production planning in an oil reservoir. Computers & Chemical Engineering, 55(0): 109-125. Teare, M., Cruickshank, R., Miller, S., Overland, S. and Marsh, R., 2013. Alberta’s Energy Reserves 2013 and Supply/Demand Outlook 2014–2023, Energy Resources Conservation Board, Government of Alberta: Calgary, Alberta. Vanegas, J.W., Cunha, L. and Deutsch, C.V., 2011. Proxy models for fast transfer of static uncertainty to reservoir performance uncertainty. AAPG Memoir(96): 203-215.
24
Table 1. Fixed DA and SP dimensions Parameter
Value (m)
DA length
850
DA width
600
SP length
150
SP width
50
DA–SP spacing
50
Table 2. First stage optimization results Varialbe
Value
θ (rad)
0.7229
∆x (m)
-64.45
∆y (m)
-273.1
Available Bitumen (NCB)
16720
Table 3 First stage optimization results using decomposition scheme 1 Optimization results
Area 1
Area 2
Area 3
θ (rad)
0.6283
0.7362
0.5642
∆x (m)
71.05
25.25
-66.45
∆y (m)
386.2
-169.5
-395.6
Available Bitumen
4056
5646
5970
Total Available Bitumen
15681
25
Table 4. First stage optimization results using decomposition scheme 2 Optimization results
Area 1
Area 2
Area 3
θ (rad)
-0.9675
0.9362
0.6429
∆x (m)
103.55
108.75
-178.45
∆y (m)
-421
226.7
59.4
Available Bitumen
4082
7063
5802
Total Available Bitumen
16947
Table 5. Athabasca reservoir characteristics (Guo et al., 2013) Parameter
Value
Parameter
Value
w (m)
75
a
1.9
k (μm2)
2.5
b
0.66
α (m2 s-1)
7.06×10-7
β
0.264
m
4
γ
0.132
νs (mm2 s-1)
10
ur
0.556
S0
0.8
ud
0.864
ϕ
0.35
vr
2.9
h (m)
26.5
vd
2.125
Table 6. Reservoir characteristics Parameter
Value
DAcell
510
dp (m)
90×10-6
Av (m-1)
6.67×10+4
τ (%)
81
26
Table 7. Design and economic parameters (Giacchetta et al., 2015; JacobsConsultancy, 2012; Lawal, 2014) Parameter
Value
Parameter
Value
DAwell
6
Discount Rate (%)
6
L (m)
850
Tax Rate (%)
30
SOR
2.7
OP (CAN$ m-3)
719.13
3
-1
-3
Basic CPF Capacity (Mm yr )
3.5
FC (M$ m )
3.023×10-7
SAGDlifetime (yr)
25
CCCPF (M$)
2480
Royalty (%)
10
CCDA (M$)
78
Table 8. Solution report for each case study c2
CPU
Optimality
(%)
time (s)
gap (%)
50
3600
1.56
100
804
0.98
150
0.28
0.78
200
0.14
0.03
250
0.02
0.01
300
0.09
0.45
27
Figure 1. SAGD well illustration (left: 3-D view of a SAGD well pair; right; 3-D view of a drainage area with 3 well pairs)
Figure 2. Overview of Optimization Framework 28
Figure 3. Net Continuous Bitumen Calculation
Figure 4. Total NCB Calculation of a DA from aerial view
29
Figure 5. Leased area with surface restriction of the reservoir for well placement
Figure 6. Porosity distribution of the leased area
30
Figure 7. Permeability distribution of the leased area
Figure 8. Net Continuous Bitumen profile
31
Figure 9. Compact DA placement plan
Figure 10. Compact DA placement plan using decomposed area (scheme 1)
32
Figure 11. Compact DA placement plan using decomposed area (scheme 2)
Figure 12. Predicted oil production rates for the Athabasca reservoir
33
Figure 13. Case study 50%–solution
Figure 14. Case study 150%–solution 34
Figure 15. Case study 200%–solution
Figure 16. Optimal NPV at different scenarios for CPF capacity 35
Figure 17. Planning developments on the compact DA placement plan
Figure 18. Oil production rate of DAs over lifetime of the project
36
Highlights
Strategic planning of oil sands development through SAGD process Two-stage optimization framework to plan the development of drainage areas Maximization of the recoverable bitumen through drainage area layout optimization Maximization of the net present value through mixed integer linear optimization
37