© H. Liu et al., published by IFP Energies nouvelles, 2020

1 Introduction

Reservoir simulators have been developed to simulate real production processes for oil, water, and gas. A simulator models the flow behavior of these fluids in porous media. It also models injection and production wells. Simulators are important tools to the petroleum industry, and have been applied to oil and gas production for decades. A simulator is applied to design a production plan before the plan is implemented in a petroleum reservoir. It is also employed to optimize oil and gas production.

Reservoir simulation is an interdisciplinary research subject that has been studied for decades, where petroleum engineering, applied mathematics, numerical methods, and computer sciences are heavily applied. There are three major types of reservoir models: black oil, compositional, and thermal. The standard black oil model considers water, oil and gas where oil stays in the oil phase, water stays in the water phase and gas can be in the oil and gas phases. The standard compositional model has one water and several hydrocarbon (oil and gas) components. This model is much more complex than the black oil model. The standard black oil and compositional models do not consider a temperature change. The thermal model is the same as the compositional model except that a temperature change is considered to model in-situ combustion and steam injection. Many researchers have been working on reservoir models, correlations for various properties, numerical methods and their efficient implementations [1–7]. Each reservoir model has its own mathematical description and numerical methods. When more features are considered, a model may be more difficult to compute. Efficient numerical treatments should be designed, such as upstream, decoupling [8] and special preconditioners like Constrained Pressure Residual (CPR) [1, 6]. When a reservoir model is large, the simulation time can be very long, such as days or even weeks. Parallel computers are employed to accelerate the reservoir simulation. A well designed parallel simulator can accelerate reservoir simulations thousands of times faster than serial simulators [1, 5].

In reservoir simulation, well modelling is complex. Each well type has its own special treatment. For example, when an injection well is modelled, the total mobility should be applied, and when a production well is modelled, only a phase mobility is required. Each well operation (constraint) has one equation, and the equations are different, such as a pressure equation for bottom hole pressure operation, rate equations for rate operations, and heating equations (heat rate and temperature) for different heating wells. In real production processes, many fluid and well operations can be applied to these wells. Various well modelling methods have been proposed [9–12]. Van Poolen et al. made early contributions to well modelling [12, 13]. Peaceman proposed most widely used well equations for square grid blocks and non-square grid blocks [9, 14, 15], which established relationships among a well rate, grid block pressure, bottom-hole pressure, a well equivalent radius, a grid size, a skin factor, permeability, and a well index. In [9], the proposed method was based on the finite difference method on square grid blocks, which is the first comprehensive study for single phase on square grid blocks [13]. Later, improved well modelling was developed for non-square grid blocks and anisotropic media [14]. Chen et al. reviewed computations of well indices and equivalent radii for horizontal and vertical wells [13, 16]. Abou-Kassem and Aziz proposed analytical methods for calculations of a well equivalent radius for an off-center well location [17]. Nolen introduced well treatments and extensions [18, 19]. A well may have multiple perforations, and each perforation is modelled individually. An overall well rate is the sum of perforations [13]. These methods can be applied to thermal reservoir simulation directly. Enthalpy is required to be calculated with well rates [20]. Isenthalpic calculations are mandatory for injection wells, which determines the status of injected fluids, such as density, pressure, phase status and composition [13, 20]. Multi-segment well models have been studied by researchers, which can treat wells with complex geometry and advanced techniques, such as tubing. Holmes et al. proposed a multi-segment well model for black oil equations [21], where each segment had its own variables, such as pressure and saturation [19]. Jiang studied multi-segment wells systematically and implemented his method in GPRS [19]. Dong applied a multi-segment well model to thermal reservoir simulation [20] to handle complex wells and advanced well treatments. Xiong developed a standalone wellbore simulator for thermal reservoir models [22]. CMG STARS manual is a great reference for thermal well modelling, which introduces various methods and extensions [10]. It also represents how the industry uses the well modelling.

This paper will review well modelling for the thermal reservoir model, including well types, well operations, well indexes and special numerical treatments. The well modelling methods can be applied to other reservoir models directly, such as the black oil and compositional models. Numerical results are presented to validate these methods and their implementations by comparing with the commercial simulator Computer Modelling Group (CMG) STARS. The structure of the paper is as follows: In Section 2, the thermal reservoir model is briefly introduced. In Section 3, well types, including injection wells, production wells and heaters, and well operations for the thermal model, including pressure control, rate control and heating control, are presented, and their mathematical models are provided. In Section 4, numerical experiments are carried out to study these well operations. The results are compared with those from the commercial simulator CMG STARS. Furthermore, scalability tests are presented to demonstrate our parallel thermal simulator and well modelling.

2 Mathematical model

Many thermal reservoir simulators have been developed by service companies and academia, and most of them have similar mathematical models [10, 20, 23, 24]. Here the main equations are introduced, which are in the standard textbooks [10, 20, 23, 24].

2.1 Mass conservation equations

In the black oil model, mass is applied to mass conservation equations. However, mole fractions are widely used by the compositional and thermal models, which are multi-phase multi-component systems. The difference is that the temperature in the compositional model is constant while the temperature changes in the thermal model. The mole fraction, x _c,α , of a component is the ratio of its mole number to the total mole number of the multi-component system in phase α (water, oil, and gas),(1)where n _c,α and n _α are the mole number of a component and the total mole number in α phase, respectively. In the thermal reservoir model, the water, oil and gas phases are considered. A component may exist in any phase. For example, the water component can exist in the water and gas phases, heavy oil stays in the oil phase only, light oil exists in the oil and gas phases, and a non-condensable gas component stays in the gas phase only. The mass conservation law of component c by its mole number is written as [25],(2)where ρ _α is the molar density of phase α, S _α is saturation, x _c,α is a molar fraction, and q _α,well is a well rate. The well rate can be negative or positive depending on the well types. The equation defines a fluid velocity among pressure, gravity, viscosity, permeability and relative permeability. For the water, oil, and gas phases, Darcy’s laws are described as,(3)where k _rα is the relative permeability of phase α, μ _α is the viscosity of phase α, p _α is the phase pressure, γ _α is gravity and z is depth. As for the relative permeability of oil k _ro, several models are available [26–29]. Here Stone’s model II [30] is applied.

2.2 Energy conservation equation

A reservoir can have many wells. Some of them inject a fluid to the reservoir, such as steam, water, solvent and gas, while others heat the reservoir. Production wells produce oil, gas and water. The wells and other factors change the energy of a reservoir, especially when hot steam is injected. The energy conservation law is introduced to model an energy change in the thermal model [10, 20, 23–25],(4)where U _α is the volumetric internal energy for phase α, U _r is rock energy capacity, H _α is enthalpy of phase α, K _T is heat conductivity of the reservoir and fluid, and Q _loss is heat loss to overburden and underburden. A simple linear correlation is applied [10] to model heat conductivity,(5)where K _w, K _o, K _g, and K _r are thermal conductivities for the water phase, oil phase, gas phase, and rock. In some models, the solid phase is also considered, and K _s is introduced for the solid phase. The heat loss is modelled by a semi-analytical method [31], which was proposed by Vinsome and Westerveld. In real application, a reservoir model can have several rock types. In this case, different parameters are used to different rock types, such as U _r.

2.3 Density

The gas phase density is modelled as,(6)where Z is compressibility factor. The factor is calculated by the RK EOS [10].

In the thermal model, the water and oil components in a liquid phase use the same density formula, which is a function of pressure, p, and temperature, T,(7) (8)where ρ _ref is the reference density at the reference temperature T _ref and reference pressure p _ref, and c _p , c _t1, c _t2 and c _pt are input parameters. The oil phase mole density, ρ _o (p, T, x _i ), is calculated by the ideal mixing rule,(9)

2.4 Viscosity

The viscosity of an oil component can be very high when temperature is low, especially heavy oil, whose viscosity can be 10⁶ cp at the standard conditions. A liquid component, such as water and oil components, has a widely accepted correlation for computing viscosity [10, 13],(10)where μ _a and b are input parameters. They are also written as avisc (μ _a ) and bvisc (b).

Since the water phase has water only, its phase viscosity can be computed right away. However, the water phase can have more than one component when salt and polymer are considered. The viscosity of the oil phase is computed by the following logarithmic mixing rule:(11)which is the same as,(12)where i is the i-th component.

The calculation of the gas component viscosity is different. Also, the gas component viscosity is not very sensitive to temperature. Another correlation is applied,(13)where μ _a and b are input parameters. They are also written as avg (μ _a ) and bvg (b). The gas phase viscosity is computed as in CMG STARS [10],(14)where M _i is the molecular weight of the i-th component.

2.5 Energy

Enthalpy is a measurement of energy in a multi-component multi-phase reservoir, which includes fluid energy, rock energy, heat loss and energy introduced by a well. Several methods have been proposed to model enthalpy, such as gas-based, liquid based, and simple Hvap [10]. The following correlations are the same as in CMG STARS.

The gas-based method is introduced here. The enthalpy of a gas component in the gas phase is computed as follows [10]:(15)where cpg1, cpg2, cpg3, cpg4, and cpg5 are coefficients for the gas component.

For the liquid phase (oil and water phases), water, heavy and light oil exist in the liquid phase. Vaporization is considered, which models an energy change when a liquid component moves from the liquid phase to the gas phase. The following formula is applied:(16)where T _crit is critical temperature of the component, and hvr and ev are input parameters, which can be measured from the lab. hvr is a non-negative coefficient. For heavy oil components, since they exist in the oil phase only, hvr is zero. The enthalpy of a liquid component in the liquid phase equals,(17)

The enthalpy of a phase is computed by a linear mixing rule. For example, the oil phase enthalpy is calculated as,(18)where H _o,i is the liquid enthalpy of the i-th component in the oil phase, and the enthalpy of the gas phase is calculated as,(19)where H _g,i is the enthalpy of the i-th gas component in the gas phase. The gas component can be water, light oil or non-condensable gas components.

The internal energy of a phase is defined as [10],(20)and the internal energy of rock is defined as,(21)whose units are volumetric units, such as J/m³. Different correlations were also introduced [13].

3 Well modelling

An injection well and a production well can have lots of perforations, which allow a fluid to enter or to be produced. Each perforation has its own properties, such as water rate, oil rate, gas rate, well index, and pressure difference control. These properties are either provided by model input or calculated by internal formulas. Peaceman’s model is applied to model well performance. For a well perforation m, its fluid rate for phase α, Q _m,α , is calculated as [9],(22)where W _I is the well index, ρ _α is fluid density at wellbore, k _rα is relative permeability, μ _α is viscosity, p _b is well bottom hole pressure defined at the depth of z _bh, p _α is the phase pressure of a grid block that contains the perforation, z is the depth of the grid block, g is the gravity constant, and γ _α is density.

In simulation implementation, is defined together, which is noted as mobility. The mobility of injection well uses total mobility,(23)which have been applied by many simulators [10, 20]. The mobility of production well is the regular mobility, , whose value is from the grid block that contains the perforation. The relative permeability and viscosity is either from last time step, which is called explicit mobility weight, or from last Newton iteration, which is called implicit mobility weight.

The treatment of density, ρ _α , is tricky. For production well, the density is from the grid block, which is straightforward. For injection well, the density is from wellbore, where flash calculation is required to obtain accurate pressure, temperature, phase distribution of injected fluid, and density at the perforation. Isentropic flash calculation is required, which treats the enthalpy of injected fluid as constant and compute pressure, temperature, phase distribution and average density.

The explicit and implicit mobility weights have been mentioned above. Another method is to utilize unweighted mobility weight, where well rate is defined as,(24)where W _I is user input value. This weight is only valid for injection well.

3.1 Well index

Well index relates well rates and grid properties, which couples wells and reservoir. Its definition is very complex, especially when the grid is complex. Well index can be set by a user reservoir model, which provides a well index for a perforation, or be calculated using analytical methods.

For analytical methods will be introduced: GEO, GEOA, KH and KHA. Here we assume the grid is structured, the grid block is hexahedron. The definitions for well index in x, y, and z directions are similar, so only z direction is introduced. Also, well index for non-parallel direction is computed similarly, whose details can be read from CMG manual [10].

A vertical well is parallel to the z direction, as shown in Figure 1. The figure shows that the well perforates a grid block, and the lengths of the grid block along the x, y, and z directions are h _x , h _y and h _z , respectively. In a reservoir, the permeability of the grid block along the x, y, and z directions are k _x , k _y and k _z , respectively. The well radius is r _w . Its well index for GEO and GEOA models is defined as,(25)where r _e is an equivalent radius, s is a skin factor, ff is well partial completion (dimensionless), frac is fraction of a circle (dimensionless). For KH and KHA models, the well index is defined as,(26)where k _h is either from user input or computed from internal formula.

For GEOA and KHA models, the equivalent radius, r _e , is defined as,(27)where h _x and h _y are the lengths of the grid block in x and y directions. For GEO and KH models, the equivalent radius, r _e , is defined as,(28)where geofac is a user input.

3.2 Pressure difference

To compute the well rate, the pressure difference, δp, between grid block and wellbore should be obtained,(29)

As mentioned above, the well bottom hole pressure, p _b , is defined at a reference depth, z _bh, and the grid block pressure is from the grid block. After each Newton iteration or time step, the grid block pressure and well bottom hole pressure are computed. Also, well bottom hole pressure is always an unknown to be computed by Newton method.

To get the pressure difference, the pressure at the wellbore should be obtained, which is equal to (p _b  – γ _α g(z _bh − z)). There are two methods. The first one is to calculate the gravity difference between two neighbor perforations using average density, where the average density may be different perforation by perforation. The second method is to use one average density for all perforation, which is easy and is friendly to parallel computing. The average density for the second method can be obtained by well rate average,(30)where m is the perforation set of a well, ρ _m,α is fluid density for phase α(w, o, g) at the m-th perforation, and Q _m,α is well rate for phase α(w, o, g) at the m-th perforation. Then the pressure difference is calculated as,(31)

3.3 Well constraints

Well operations are complex. A well may switch its operation dynamically, such as pressure control, rate control and temperature control. Its type may be changed too. For example, an injection well may serve as production well at certain time. This section presents common used well constraints (operations) with their mathematical descriptions.

3.3.1 Fixed bottom hole pressure

Fixed bottom hole pressure well operation means to apply constant wellbore pressure at a given reference depth, whose mathematical model is written as,(32)where c is a constant.

3.3.2 Fixed rate

Fixed rate well constraint means to apply constant rate to a well. For example, if an injection well injects 100 bbl water per day to a reservoir, its well constraint is fixed water rate. Other rate controls are also applied to well operation, such as oil rate, water rate, gas rate, liquid rate (oil and water) and fluid rate (water, oil and gas). Also the rate can be measured at reservoir condition or surface condition. When surface condition is applied, flash calculation is required to calculate oil, water and gas phases distribution at surface condition, which uses segregated flash calculation or PT-flash calculation.

For a specific phase α, such as water, oil or gas, its fixed rate constraint is written as,(33)where c is a constant value. Fixed liquid rate constraint is written as,(34)

Fixed fluid rate is written as,(35)

Here are some widely applied well rate constraints for thermal reservoir simulations: (1) STW: water rate, surface condition; (2) STO: oil rate, surface condition; (3) STG: gas rate at surface condition; (4) STL: liquid rate, surface condition; (5) STF: total fluid rate, surface condition; (6) BHW: water rate, reservoir condition; (7) BHO: oil rate, reservoir condition; (8) BHG: gas rate, reservoir condition; (9) BHL: liquid rate, reservoir condition; (10) BHF: total fluid rate, reservoir condition; (11) STEAM: live steam rate, reservoir condition, cold water equivalent.

3.3.3 Constant heat transfer model

There are some different heater models, such as constant heat transfer model, convective heat transfer model and heat well. The first two types can be defined in any grid block. However, the heat well (HTWELL in CMG STARS) can only be defined in a real well, such as injection well and production well. The constant heat transfer model means in some grid blocks, there exists heat transfer at constant rate, such as 1000 Btu/day,(36)

3.3.4 Convective heat transfer model

Constant heat transfer model simulates constant heat exchange while convective heat transfer model defines dynamic heat transfer, which is controlled by two parameters: uhtr (proportional heat transfer coefficient) and tmpset (temperature setpoint). If uhtr is positive, it means to gain heat from source, the heat rate in a grid block is calculated as,(37)

If uhtr is negative, it means the reservoir loses heat, the heat rate in a grid block is calculated as,(38)where T is the reservoir temperature.

3.3.5 Heater well

As mentioned above, the constant and convective heat transfer models can be defined in any grid block. Another heat model is also developed in CMG STARS and our simulator, which is noted as HTWELL as in CMG STARS. This type of heater is defined in a production or injection well, which has the same perforations as the well contains the heater well. This heater well is more complicated than constant and convective heater transfer models, which has more controls, such as heat rate model (HTWRATE or HTWRATEPL in CMG STARS), temperature model (HTWTEMP in CMG STARS), heat index model (HTWI in CMG STARS), and dual rate/temperature model. The dual rate/temperature model has two direction controls: uni-directed (UNIDIRECT in CMG STARS) and bi-directed (BIDIRECT in CMG STARS).

For heat rate control (model), the heat rate in a perforation m is calculated as,(39)where q is the heat rate, Q _hspec is total heat rate defined by HTWRATE, L _m is the length of the layer well completion, and L _w is the total well length (sum of L _m ).

For temperature model, the heat rate in a perforation is calculated as,(40)where I _m is the heat conduct index (or heat index), T _wspec is specify wellbore temperature, T _m is grid block temperature. We should mention that there are two methods for calculating heat conduct index: (1) use thermal conductivity formula introduced in mathematical model section; (2) use heat index introduced here (by turning HTWI on in CMG STARS). For heat index model, if HTWI is turned on, the user input well index or internal index can be converted to heat index. Also, in our simulator, heat conduct index for HTWELL allows finer definitions.

When dual rate/temperature model is enabled, the rate model and temperature model are switched automatically. For heating (Q _hspec > 0), the heat rate in a layer is defined as,(41)where ΔT _m is defined as,(42)

The T _k is reservoir temperature in a grid block. The UNIDIRECT option shuts down heater when temperature difference is zero; while BIDIRECT allows heating and cooling (heat loss).

For cooling (Q _hspec < 0), the heat rate in a layer is defined as,(43)where ΔT _m is defined as,(44)

The UNIDIRECT option shuts down cooling well, and the BIDIRECT option allows bidirectional heat transfer. In both cases, the BIDIRECT can simulate autoheater and autocooler.

3.3.6 Autoheater and autocooler

Autoheater and autocooler work with constant heater and convective heater. When autoheater is applied to a location (UBA, User Bock Address), if reservoir temperature is below a temperature set (tmpset), the heat gain is from the minimum of constant heater and convective heater.

Autocooler is applied as automatic cooler to mimic heat loss. If autocooler is turned on, the heat loss is the minimum of constant heater and convective heater. If reservoir temperature is above a temperature set, heat is being lost. If it is equal to or below a temperature set, constant heater and convective heater are turned off. Autocooler has the same input as autoheater. However, if both autoheater and autocooler are applied to a grid location, autocooler is ignored.

3.3.7 Subcool control

Subcool control is also known as steamtrap, which is used to prevent the production of live steam. It does this by keeping the well’s flowing bottomhole pressure (and hence the pressure in the grid block containing the well) high enough that live steam does not appear in the well block [10].

The well constraint equation solved is written as,(45)where c is a pre-defined temperature difference, T _sat is the steam saturation temperature corresponding to wellbore pressure p _wb defined in the perforation, and T _k is the temperature defined in the grid block that contains the perforation.

4 Numerical studies

This section presents the numerical results, which has two sections. One section studies validation of different well constraints, which are compared with CMG STARS. The second section gives scalability results.

The benchmark models based on some standard models, such as SPE4 [32]. The parameters are modified to test those different types of well constraints, reservoir properties, and different chemicals (heavy oil, light oil and non-condensable gas). Equivalent models for CMG STARS are created and numerical results are compared.

The SPE4 model is a 3-D model with a mesh grid 9 × 9 × 4. The grid sizes of x-direction and y-direction are both 29.17 ft. The grid size of z-direction is 10 ft. The porosity is a linear function of the pressure. It has a constant porosity 0.3 at the reference pressure 75 psi. The heat capacity of the rock is 35 Btu/(ft³ F). It has constant permeabilities of 1000, 1000, 500 md in x, y and z directions, respectively. The thermal conductivities of the rock, water, oil and gas are all 24 Btu/ft day c F. The reference pressure of the reservoir and the surface are both 14.7 psi. The reference temperature of the reservoir and the surface are both 60 F. The initial water mole density is 3.464. The compressibility of the water is 3.999e-6 1/psi. The mole weight of the water component is 18.02 lb/lbmol. The critical pressure and critical temperature of the water are 3206.2 psi and 705.4 F respectively. The volumetric heat capacities of formation adjacent to the reservoir at the overburden and underburden are 35 Btu/ft³ × F. The thermal conductivities of formation adjacent to the reservoir at the overburden and underburden are 24 Btu/ft day F.

4.1 Well constraints

There are several well constraints implemented in our simulator such as the fixed bottom hole pressure, constant heat transfer model, convective heat transfer model, subcool control and heater constraints. Here we compare our results of the SPE4 model with different well constraints to the results from CMG STARS.

4.1.1 Fixed bottom hole pressure

To test the fixed bottom hole pressure constraint, the bottom-hole pressures of the production wells are fixed to be a constant.

Example 4.1

There are heavy oil and light oil in the reservoir. The oil mole densities of the heavy and light oil are 0.10113 and 0.2092 lb/lbmole, respectively. The compressibilities of both the heavy oil and light oil are 5.0e-6. The mole weights of the heavy oil, light oil are 600, 250 lb/lbmol. The critical pressure and critical temperature of the heavy oil are 0 psi and 0 F, respectively. The critical pressure and critical temperature of the light oil are 225 psi and 800 F respectively. The volume of the rock is assumed constant which preserves the cell rock mass and hence rock heat capacity at a given T. The initial condition is as follows: The water and oil saturations are 0.45 and 0.55, respectively. The mole fractions of the heavy oil and the light oil are 0.6 and 0.4, respectively. The initial temperature is 125 F. The initial pressure is 4000 psi. The overburden temperature is 125 F. There are one injection well and four production wells P _j (j = 1,⋯, 4) in the reservoir. All the wells are vertical with perforations from the first layer to the 4th layer. The injection well locates at (5, 5) with well index 2.0e4 md ft. The production wells locate at (1, 1), (1, 9), (9, 1) and (9, 9) grid positions with well indices 2.0e4, 3.0e4, 4.0e4, 5.0e4 md ft. The well constraints are set as follows: the maximum water injection rate is 100 bbl/day at the temperature 450 F with 30% steam; the minimal bottom hole pressures of all the four production wells are set to 17 psi. The simulation time is 730 days.

The bottom hole pressures of the injection well and the production well P ₁ are shown in Figure 2. From this figure, we can see that our bottom-hole pressures of the injection well and well P ₁ match the CMG STARS results exactly. The oil rate, the gas rate and the water rate are also shown in Figure 3. The oil rate is plotted in Figure 3a from which we can see that our result almost overlaps with the result from CMG STARS. The gas rate and the water rate are depicted in Figures 3b and 3c which show that during the first 200 days, both the gas rate and the water rate from our simulator are same as the CMG STARS. After 200 days, although there are a little difference, our rates are still close to the CMG STARS. The field rates of the four production wells are shown in Figure 4.

4.1.2 Constant heat transfer model

In this example, the constant heat transfers are added to the production wells to validate the implementation.

Example 4.2

All the parameters are same as Example 4.1 except that the overburden temperature equals to the initial reservoir temperature 60 F and the initial mole fraction of heavy oil is 1.0. The well constraints are kept same as Example 4.1 except that for all the production wells, the constant heat transfer rates are added such that heat exchanges at some perforations. For every production well, the heat transfer rate is 1.0e5 Btu/day at every well perforation.

The bottom-hole pressure and the injection rate of the injection well are given in Figure 5. From this figure, we can see that the result of the bottom-hole pressure of the injection well from our simulator is exactly same as the result from CMG STARS. Although the injection rate oscillates, the relative error is lower than 0.01%. The oil production rates and the water production rates of well P ₁ are shown in Figure 6. The field oil production rates and the field water production rates of all the production wells are depicted in Figure 7. These figures show that the results of our simulator match the results from CMG STARS very well.

4.1.3 Convective heat transfer model

In this example, the convective heat transfers are added to the production wells. The numerical results from our simulator are compared to the results of CMG STARS for validation.

Example 4.3

There is only heavy oil in the reservoir. The oil mole density of the heavy is 0.10113 lb/lbmole. The compressibility is 5.0e-6. The mole weight is 600 lb/lbmol. The critical pressure and critical temperature of the heavy oil are 0 psi and 0 F respectively. The volume of the rock is assumed constant which preserves the cell rock mass and hence rock heat capacity at a given T. The initial condition is as follows: The water and oil saturations are 0.45 and 0.55, respectively. The mole fraction of the heavy oil is 1.0. The initial temperature is 125 F. The initial pressure is 75 psi. The overburden temperature is 125 F. There are one injection well and 2 production wells P _j (j = 1,⋯, 2) in the reservoir. All the wells are vertical with perforations from the first layer to the 4th layer. The injection well locates at (1, 1) with well index 2.0e4 md ft. The production wells locate at (9, 1), (5, 5) grid positions with same well index 1.0e4 md ft. The well constraints are as follows: the maximum water injection rate is 50 bbl/day at the temperature 450 F with 0% steam; the minimal bottom hole pressures of the two production wells are set to 17 psi. For all the three wells, the convective heat transfer rates are added such that the wells can gain heat from the sources. For the injection well, the gain coefficient of a temperature controller is 1.0e6 Btu/day F at every well perforation. For the two production wells, the gain coefficient of a temperature controller is 4.0e4 Btu/day F at every well perforation. The temperature setpoints are all equal to 500 F. The simulation time is from the first day to the 730th day.

Similar to Example 4.2, first the bottom-hole pressure and the injection rate of the injection well are shown in Figure 8. From the left sub-figure we can see that the bottom-hole pressures from our simulator and from CMG STARS overlap. The right sub-figure shows that the injection rates of the injection well from both our simulator and CMG STARS oscillate. However, amplitudes of the oscillations are both less than 0.05%. Hence, the results are still very close. The oil production rate and the water production rate are depicted in Figure 9 which shows the good match of our results to the results of CMG STARS. The field production oil rate and the field production water rate are given in Figure 10 from which we can see that our results are exactly same as the results from CMG STARS.

4.1.4 Subcool control

Example 4.4

There is only heavy oil in the reservoir. The oil mole density of the heavy is 0.10113 lb/lbmole. The compressibility is 5.0e-6. The mole weight is 600 lb/lbmol. The critical pressure and critical temperature of the heavy oil are 0 psi and 0 F, respectively. The volume of the rock is assumed constant which preserves the cell rock mass and hence rock heat capacity at a given T. The initial condition is as follows: The water and oil saturations are 0.45 and 0.55, respectively. The mole fraction of the heavy oil is 1.0. The initial temperature is 125 F. The initial pressure is 75 psi. The overburden temperature is 125 F. There are one injection well and 2 production wells P _j (j = 1,⋯, 2) in the reservoir. All the wells are vertical with perforations from the first layer to the 4th layer. The injection well locates at (1, 1) with well index 1.0e4 md ft. The production wells locate at (9, 1), (5, 5) grid positions with same well index 1.0e4 md ft. The well constraints are as follows: For the injection well, the maximum injection water rate is 30 bbl/day at the temperature 450 F with 60% steam and the maximum bottom-hole pressure 1500 psi. For the two production wells, the minimum bottom-hole pressures are both 17 psi with steam trap value 20 F. The simulation time is from the first day to the 2190th day.

The bottom-hole pressure and the injection rate of the injection well are plotted in Figure 11. From the left sub-figure, we can see that the bottom-hole pressure of the injection well from our simulator is exactly same as the one from CMG STARS. From the right sub-figure, we can see that the injection rates from our simulator and from CMG STARS are both very close to 30 bbl/day although both of them have some small oscillations. For the production well P ₁, the oil production rate and the water production rate are report in Figure 12 from which we can see that there is little difference between the results from our simulator and CMG STARS. For all the production wells, the field production oil rate and the field production water rate are given in Figure 13. From its left sub-figure, we can see that our result of the oil production rate are almost same as the result from CMG STARS. Same happens to the water rate shown in the right sub-figure.

4.1.5 Heater constraints

Example 4.5

There is only heavy oil in the reservoir. The oil mole density of the heavy is 0.10113 lb/lbmole. The compressibility is 5.0e-6. The mole weight is 600 lb/lbmol. The critical pressure and critical temperature of the heavy oil are 0 psi and 0 F, respectively. The volume of the rock is assumed constant which preserves the cell rock mass and hence rock heat capacity at a given T. The initial condition is as follows: The water and oil saturations are 0.45 and 0.55, respectively. The mole fraction of the heavy oil is 1.0. The initial temperature is 125 F. The initial pressure is 75 psi. The overburden temperature is 125 F. There are one injection well and two production wells P _j (j = 1,⋯, 2) in the reservoir. All the wells are vertical with perforations from the first layer to the 4th layer. The injection well locates at (1, 1) with well index 1.0e4 md ft. The production wells locate at (9, 1), (5, 5) grid positions with same well index 1.0e4 md ft. The well constraints are as follows: For the injection well, the maximum bottom-hole pressure is 5000 psi with the maximum injection water rate 300 bbl/day. Both of the production wells are heater wells. For well P ₁ , the minimum bottom-hole pressure is 17 psi with maximum heating rate 1.0e6 Btu/day and wellbore temperature 600 F. The temperature model is used. The thermal conductivity is computed internal instead of using the well index of the perforation. The heat transfer for the well P ₁ at every perforation is 1.0e6 Btu/day. The gain coefficient of a temperature controller is 4.0e4 Btu/day F at every well perforation. The temperature setpoint of the temperature controller is 500 F. For the well P ₂ , the minimum bottom-hole pressure is 17 psi with maximum heating rate 3.4e6 Btu/day and wellbore temperature 611 F. The dual rate/temperature mode is used. The thermal conductivity is computed internal instead of using the well index of the perforation. The heat transfer for the well P ₂ at every perforation is 1.0e6 Btu/day. The gain coefficient of a temperature controller is 4.0e4 Btu/day F at every well perforation. The temperature setpoint of the temperature controller is 500 F. The simulation time is from the first day to the 730th day.

Similarly, first the bottom-hole pressure and the injection rate of the injection well are given in Figure 14 which shows the results of both the bottom-hole pressure and the injection rate are almost the same as the CMG STARS. The oil production rate and the water production rates of well P ₁ plotted in Figure 15 have a little difference from the results of CMG STARS. This also happens to the field production oil rate and the field water rate depicted in Figure 16.

4.2 Scalability

Two supercomputers are employed, Cedar and Archer. Cedar is a supercomputer from Compute Canada, which has a total of 58 416 CPU cores for computation, and 584 GPUs [33]. It uses Intel OmniPath network with bandwidth of 100 Gbit/s. Archer is from UK National Supercomputing Service, which is a Cray XC30 supercomputer. It has 4544 standard memory nodes, each of which has two 2.7 GHz, 12-core E5-2697 v2 (Ivy Bridge) series processors [34]. These nodes are connected by Cray Aries network in a Dragonfly topology.

Example 4.6

This example presents a large thermal model with a grid dimension of 360 × 2000 × 1600, which has 1.2 billion grid blocks. 120 computation nodes are employed using the Cedar supercomputer, and up to 960 CPU cores are used. Each node runs 2, 4, and 8 CPU cores (MPIs). The standard Newton method is applied with a tolerance of 1e-10 and maximal iterations of 10. The linear system has 6 billion unknowns, and the linear solver BICGSTAB is applied, which uses a tolerance of 1e-10 and maximal iterations of 100. The preconditioner is the CPR-FPF method [ 1 ]. This model has several wells. Some inject steam to reservoir and some produce water, oil and gas. Table 1 presents running time and memory used. Figure 17 shows the scalability.

Table 1 presents running times for simulation and linear solver, speedup and memory required. We can see that when more CPU cores are employed, the simulation time is cut, which means parallel computing is a great tool for reservoir simulations. The linear solver uses around 50% of the total simulation time, and it is clear that efficient linear solver methods are critical to the success of reservoir simulations. Also, huge amounts of memory (around 10 000 GB) are required for large-scale reservoir simulations, which is impossible for our desktop computers. Figure 17 shows the simulator, linear solver and preconditioner have excellent scalability. This example indicates our parallel simulator can handle large-scale models, and the linear solver and preconditioner can solve linear systems with billions of unknowns.

Example 4.7

This example studies a larger thermal model with a grid dimension of 9.2 billion grid blocks and the resulted linear systems have 46 billion unknowns. 2048 nodes are employed. The Newton method is applied with a tolerance of 1e-4 and maximal iterations of 10. The linear solver is BICGSTAB with a tolerance of 1e-3 and maximal iterations of 100. The preconditioner is the CPR-FPF method with GJE decoupling. Table 2 presents running time and memory used, and Figure 18 shows the scalability (overall speedup) curve.

Table 2 shows summary of this example. In this example, each node uses 2, 4, 6 and 12 cores (MPIs). When multiple cores are employed on one node, these cores compete memory bandwidth and their effective memory bandwidth is reduced. Even though, from the table, we can see that when more cores are employed, the running time is cut, and good acceleration is observed. When 4 and 6 cores are employed, the linear solver and simulator have efficiency around 90%, which is excellent for parallel computing. When 12 cores are employed, the efficiency for solver is still over 80% and the overall efficiency is 76%, though effective bandwidth is reduced. Around 90 000 GB are required for this billion size thermal simulation. Figure 18 compares ideal scalability and actual scalability. We can conclude that our parallel thermal simulator can compute extremely large-scale thermal models and it has excellent scalability. It can handle well modelling in large-scale reservoir models.

Example 4.8

This example studies a large thermal model with a grid dimension of 20 billion grid blocks using 4096 nodes and the resulted linear systems have 100 billion unknowns. The Newton method is applied with a tolerance of 1e-4 and maximal iterations of 10. The linear solver is BICGSTAB with a tolerance of 1e-3 and maximal iterations of 100. The preconditioner is the CPR-FPF method with GJE decoupling. Table 3 presents running time and memory used, and Figure 19 is the scalability.

The overall scalability is good but the linear solver has better scalability. The memory consumption is proportional to grid blocks, which means if more computation resource is available, larger models can be studied. Since the thermal reservoir simulator has good scalability, the same model can be run faster if using more computation nodes. These examples indicate that our thermal simulator can handle giant thermal models and well modelling in large-scale parallel reservoir simulations. Meantime, we have observed that when more cores of a processor are employed, the scalability tends to reduce, which could be caused by memory bandwidth, network or algorithms, and there is room to investigate and to improve its performance.

5 Conclusion

This paper reviews well modelling for thermal reservoir model. Different well types are introduced, including injection well, production well and heaters (heating well). The Peaceman well model is applied, where mobility, equivalent radius and different well index definitions are presented, especially the well index. It has many types of definitions and some extra parameters are considered, which are practical for real applications. Various well operations and their mathematical equations are presented. These are the most commonly used methods in thermal simulations. Each operation has its own equation, and when we develop simulator, it means we have to implement different calculations and different Jacobian systems. Also, the operations have to be checked and switched during simulations. Numerical results are presented to study the well modelling and they are compared with CMG STARS. The results show they match CMG STARS very well and our in-house parallel simulator has excellent scalability, which can handle large-scale reservoir models.

Acknowledgments

The support of Department of Chemical and Petroleum Engineering, University of Calgary and Reservoir Simulation Research Group is gratefully acknowledged. The research is partly supported by NSERC/AIEES/Foundation CMG, AITF iCore, IBM Thomas J. Watson Research Center, and the Frank and Sarah Meyer FCMG Collaboration Centre for Visualization and Simulation. The research is also enabled in part by support provided by WestGrid (www.westgrid.ca), SciNet (www.scinethpc.ca) and Compute Canada Calcul Canada (www.computecanada.ca).

References

All Tables

All Figures