Energy Fuels 2011, 25, 251–259 Published on Web 12/07/2010
: DOI:10.1021/ef101225g
Parameter Space Reduction and Sensitivity Analysis in Complex Thermal Subsurface Production Processes Jacob H. Bauman and Milind D. Deo* Department of Chemical Engineering, University of Utah, 50 South Central Campus Drive, Room 3290, Salt Lake City, Utah 84112, United States Received September 9, 2010. Revised Manuscript Received November 2, 2010
As conventional resources for liquid fuels in the world become scarcer and less secure, there is a need to develop other feasible resources. Oil shale is a massive resource local to the United States for potential liquid fuel production. In situ oil shale processing strategies are attractive for reduced environmental impact (in comparison to surface production operations) and provide access to resources inaccessible to mining. The efficiency of feasible and economical development is greatly enhanced with predictive power that is both efficient and accurate. However, modeling thermal subsurface processes is a complex problem, involving many simultaneously occurring physical phenomena. In this study, an oil reservoir simulator capable of representing thermal processes was used to explore the impact and interplay of various pertinent parameters to an in situ oil shale processing strategy. A statistical methodology was developed using designed factorial experiments (simulations) to expose probable dominating parameters, including synergistic or diminutive interactions between parameters. An empirical regression model or response surface was built from the simulated data. Monte Carlo simulations were used to characterize the response surface and to estimate the uncertainty in predicted oil recovery results because of the explored parameters.
In situ thermal processing is a complex process that requires understanding of multiple phenomena at multiple scales. The thermal transformation of organic matter (kerogen in oil shale) to useful fuels (liquids and gases) requires a detailed understanding of the parent composition, the transformation pathways, and the products. The organic matter coexists with inorganic and mineral matter, and the heat and mass transfer at the grain scale affect process effectiveness.
Introduction The world is facing several interesting energy challenges. Conventional liquid fuel resources are becoming scarcer. Carbon dioxide emissions associated with global warming demand attention and technological solutions. Secure energy sources to supply increasing global demand will also be necessary to address challenges moving forward. Understanding complex thermal and reactive subsurface processes will facilitate technological development, including production of oil from oil shale and oil sands, thermal treatment of underground coal, carbon dioxide sequestration, geothermal energy production, etc. Oil shale processing technology development is attractive because of the massive resources within the United States. Resource estimates in the Green River formation located in the United States range from 1.5 trillion to 1.8 trillion barrels of original oil in place from relatively rich shales exceeding 15 gallons/ton.1,2 Two major processing strategies exist for converting oil shale to oil: ex situ and in situ. Ex situ strategies include mining organic-rich shale followed by crushing and pyrolysis heating. Various ex situ pyrolysis heating strategies exist. In situ processing strategies attempt to convert the organic matter to oil underground by some form of heating and then to produce that oil in production wells. Because heat input is required in both types of processes, heating efficiency is crucial for any successful strategy.
Description of the Important Physical Phenomena and Related Properties Heat transfer through reservoir scale systems is not trivial. A wide array of heating strategies, including resistive heating wells,3 fracture injection with conductive material,4 underground rubblization followed by well heating,5 or radio-frequency heating,6 have been proposed and developed. Effective heating of the reservoir is crucial to the efficiency of any in situ thermal processing strategy. Modeling of various strategies requires fundamental understanding of the physics, including thermal conductivities and heat capacities of inhomogeneous reservoir materials, convection, phase changes, heats of reactions, and heat losses because of inefficiencies or losses to reservoir boundaries. Estimating the (3) Wellington, S. L.; Berchenko, I. E.; Rouffingnac, E. P.; Fowler, T. D.; Ryan, R. C.; Shalin, G. T.; Stegemeier, G. L.; Vinegar, H. J. U.S. Patent 6,880,633, 2005. (4) Symington, W. A.; Olgaard, D. L.; Otten, G. A.; Phillips, T. C.; Thomas, M. M.; Yeakel, J. D. ExxonMobil’s Electrofrac process for in situ oil shale conversion. Proceedings of the 26th Oil Shale Symposium; Golden, CO, Oct 17, 2006; http://www.ceri-mines.org/documents/R05bBillSymington-rev_presentation.pdf (accessed on July 21, 2009). (5) Chevron USA, Inc. Oil Shale Research, Development and Demonstration Project Plan of Operation; Chevron USA, Inc.: Houston, TX, 2006; http://www.blm.gov/pgdata/etc/medialib/blm/co/field_offices/white_river_field/oil_shale.Par.37256.File.dat/OILSHALEPLANOFOPERATIONS.pdf (accessed on June 23, 2010). (6) Kasevich, R. S.; Kolker, M.; Dwyer, A. S. U.S. Patent 4,140,179, 1979.
*To whom correspondence should be addressed. E-mail: milind.deo@ utah.edu. (1) Dyni, J. R. Geology and resources of some world oil-shale deposits. Oil Shale 2003, 20 (3), 193–252. (2) Bartis, J. T.; LaTourrette, T.; Dixon, L.; Peterson, D. J.; Cecchine, G. Oil Shale Development in the United States: Prospects and Policy Issues; RAND Corporation: Santa Monica, CA, 2005; pp 5-9 (http://www. rand.org/pubs/monographs/2005/RAND_MG414.pdf). r 2010 American Chemical Society
251
pubs.acs.org/EF
Energy Fuels 2011, 25, 251–259
: DOI:10.1021/ef101225g
Bauman and Deo
Simulator Description
significance of these various modes of heat transfer could simplify such models because physical field data are sparse and expensive. As complex organic material, such as kerogen in oil shale, is heated, it is converted into lighter oil and gas products. A variety of kinetic transformation mechanisms have been reported.7-9 Time and temperature histories of these complex organic materials can have significant implications on product distributions of literally thousands of components with their associated properties. Typically, these components are represented with relatively few pseudo-components. The complexity of the reaction network to represent these components is a major consideration.10 Kinetic parameters, such as activation energy, frequency factor, stoichiometry, etc., greatly depend upon the complexity of the component representation. For example, in one model from Braun and Burnham, the activation energy for fractions of various representations of type-II kerogen vary from 47 to 54 kcal/mol.11 The stoichiometry of reactants and products depends upon the molecular weights and elemental representations of the lumped species to conserve mass and elemental balances. For further cracking reactions of lumped components, the kinetic parameters are dependent upon molecular weights, aromaticity, thermodynamics, etc. of all species represented by the representative component. Significant variation in organic (kerogen) and inorganic material within and between resources has implications on the appropriate kinetic representation of thermal processes as well. Permeability dynamics (initial permeability and its evolution as the process unfolds) are crucial to a successful in situ oil shale strategy. Oil shale resources are typically characterized with very low initial permeability. As solid kerogen is heated and converted to liquid and gaseous products, pore space is created. A permeability increase is possible. Experimental studies have shown expansion followed by subsidence of oil shale rock as it is heated.12 Decomposition of inorganic rock at high heats would have permeability implications, including possible microfracturing of the rock, for some heating strategies. Another possible major contributor to permeability dynamics is coke plugging of permeability pathways. Relative permeability correlations are used in reservoir simulation to account for multiphase Darcy flow in permeable rock. The relative permeability models are based on experimental data or are empirically constructed, but the accuracy in relative permeability representation in complex reservoirs may be crucial. Thermal effects on phase viscosities and dynamic capillary effects with changing rock mechanics add complexity to relative permeability representation. Significant interplay of parameters within and between these physical phenomena can exist. For example, high temperatures required for oil shale pyrolysis have significant impact on organic material composition, phase and flow behavior, and possibly, geomechanics. Although the parameters are supplied to governing equations and theoretical models in a simulator, the impact that each parameter has on the final recovery of oil predicted is not easily determined or available. This information potentially would supply researchers and modelers the parameters that are of the greatest importance and, therefore, need the greatest attention for accurate prediction of such processes.
The principal balance equations needed to solve thermal reservoir problems are species and energy balance equations. RCi ¼
Np Nf X D X ð½ ðφSp Fp xp, i ÞÞ - rðð xp, i Fp νBp ÞÞ Dt p ¼ 1 p¼1 Nr X -ð sr, i R r þ Qm, i Þ r¼1
RE ¼
Np D X ð½ ðφSp Fp Up Þ þ Urock Þ - rK c rT Dt p ¼ 1
- rð
Nf X p¼1
Nr X ^ p Fp νBp Þ - ð H ΔHr R r þ Qe Þ r¼1
The flow term is typically calculated with Darcy’s law. krp ðrPp þ γp rzÞ νBp ¼ K μp STARS from the Computer Modeling Group is capable of performing four-phase multicomponent thermal reservoir simulations.13 Equilibrium calculations are K-value-based. The simulations in this study were run with STARS. In these simulations, vertical heating wells surround a producer in a seven-point pattern. The heating wells supply only heat without injecting any fluids. Only a triangular wedge with fractions of two heaters and a fraction of one producer is discretized and calculated, as shown in Figure 1.14 The initial dimensions of the wedge were 53 ft between heaters and 50 ft thickness of the reservoir. The temperature at the heating wells was raised quickly to approximately 800 °C and then controlled at about 650 °C to supply sufficient heat for adequate heat transfer throughout the reservoir but also to maintain reasonable temperatures near the heating well. Geological data from the mahogany zone in the Uinta Basin well U059 was used to estimate the richness of the layers in the reservoir. The richness of the layers varies from 12.5 to 25 wt % of hydrocarbon material in the oil shale.15 All of the hydrocarbons were assumed to be kerogen, and this kerogen was assumed to occupy 90-95% of the pore space, with the rest of the pore space being occupied with 99% gas saturation and 1% water saturation. The initial kerogen, water, and gas volumes in the pore space will vary between resources or sections of a resource. This could have important implications on the heat-transfer dynamics dependent upon the initial water mass and volume. Kerogen was specified with a constant solid density; therefore, rich layers were assigned higher porosities and lean layers were assigned lower porosities. Porosity is defined here as the total pore space occupied by kerogen, water, and gas. Fluid porosity refers to the pore space occupied by liquids and gases. As solid kerogen is converted to liquids and gases, fluid porosity increases, while total pore space (as defined in this study) remains unchanged. Green River oil shale has been characterized with low initial
(7) Allred, V. D. Chem. Eng. Prog. 1966, 62, 55. (8) Rajeshwar, K.; Nottenburg, N.; Dubow, J. J. Mater. Sci. 1979, 14, 2025. (9) Strizhakova, Y. A.; Usova, T. V. Current trends in the pyrolysis of oil shale: A review. Solid Fuel Chem. 2008, 42, 197. (10) Pepper, A. S.; Dodd, T. A. Simple kinetic models of petroleum formation. Part II: Oil-gas cracking. Mar. Pet. Geol. 1995, 12, 321. (11) Braun, R. L.; Burnham, A. K. Chemical reaction model for oil and gas generation from type I and type II kerogen. LLNL Report UCRL-ID-114143; Lawrence Livermore National Laboratory (LLNL): Livermore, CA, 1993. (12) Mattson, E. D.; Palmer, C. D.; Johnson, E.; Huang, H.; Wood, T. Permeability changes of fractured oil shale cores during retorting. Proceedings of the 29th Oil Shale Symposium; Golden, CO, Oct 20, 2009.
(13) Computer Modeling Group. STARS User Manual; Computer Modeling Group: Calgary, Alberta, Canada, 2007. (14) Bauman, J. H.; Huang, C. K.; Gani, M. R.; Deo, M. D. Modeling of the in-situ production of oil from oil shale. Oil Shale: A Solution to the Liquid Fuel Dilemma; American Chemical Society (ACS): Washington, D. C., 2010; ACS Symposium Series, Vol. 1032, pp 135-146. (15) Vanden Berg, M. D.; Dyni, J. R.; Tabet, D. E. Utah Oil Shale Database; Utah Geological Survey: Salt Lake City, UT, 2006; OFR 469.
252
Energy Fuels 2011, 25, 251–259
: DOI:10.1021/ef101225g
Bauman and Deo
Figure 1. Visualization of a simulated section with dimensions.
fluid porosity,16 although this can vary between resources. The simulations studied in this paper have relatively low initial fluid porosities and are quite dry. Very little initial water is present in these simulations. Horizontal permeability varied from 0.1 md at the bottom of the reservoir to 1 md at the top, and vertical permeability varied from 0.05 md at the bottom of the reservoir to 0.5 md at the top. Initial permeability is typically low for oil shale, but permeability models where permeability increases as kerogen is converted to fluids would make permeability dependent upon initial kerogen richness as the process unfolds. It is likely that a permeability creation/ destruction model will be necessary to model this process because of solid kerogen conversion to fluids, volume expansion of the rock, and coke plugging of pores. Increased permeability would allow fluids to flow more efficiently through the reservoir. If significant permeable pathways develop, the residence time of products is reduced, which has compositional implications because liquid oils could further transform to gases and coke. Also, with increased permeability, there is greater opportunity for convective heat transfer for improving heating efficiency. Typical permeability creation/destruction models relate permeability to fluid porosity. In this particular study, these dynamic effects were not considered. A multiple reaction scheme was used to estimate kerogen decomposition to products. All hydrocarbons were lumped into seven representative components: kerogen, heavy oil (HO), light oil (LO), gas, methane (CH4), char, and coke. The reaction scheme was adapted from a previous study17 and is similar to other kerogen decomposition models,18 although the reaction scheme in this study is relatively simple. These representative components are lumped together based on molecular weight. Many more representative components lumped according to other physical characteristics, such as density, viscosity, aromaticity, solubility, etc., could also be represented depending upon the desired complexity in the
reaction scheme. The ability to develop a more complex set of lumped components also depends upon the available data from experiments where these physical properties of interest are analyzed. However, increased complexity can greatly increase computational cost. A simplified reaction scheme assumes that the most important physical properties for predicting reservoir behavior depend upon the molecular weight of all of the components in a real system. kerogen f HO þ LO þ gas þ CH4 þ char
ð1Þ
HO f LO þ gas þ CH4 þ char
ð2Þ
LO f gas þ CH4 þ char
ð3Þ
gas f CH4 þ char
ð4Þ
char f CH4 þ gas þ coke
ð5Þ
It should also be noted that this reaction scheme does not include mineral reactions. Carbonate minerals present in oil shale resources in the Green River formation can also decompose at the high temperatures encountered. These mineral transformations would have implications for CO2 emissions and also would have a relationship with porosity and permeability dynamics. Studies for carbonate mineral decomposition have been performed to determine the mineral reactions that could be involved.19 Experimental Design Factorial experimental designs give experimenters and analysts efficient tools to understand the impact that parameters have on a response in a process. Unlike “one at a time” experiments, factorial designs allow for the researcher to estimate interactions between parameters with fewer experiments. These factorial designs allow researchers to evaluate many factors together. Experimental design methods were primarily developed for quality assurance purposes but have
(16) Tisot, P. R. Properties of Green River oil shale determined from nitrogen adsorption and desorption isotherms. J. Chem. Eng. Data 1962, 7, 405–410. (17) Braun, R. L.; Burnham, A. K. PMOD: A flexible model of oil and gas generation, cracking, and expulsion. Org. Geochem. 1992, 19, 161– 172. (18) Stainforth, J. G. Practical kinetic modeling of petroleum generation and expulsion. Mar. Pet. Geol. 2009, 26, 552–572.
(19) Campbell, J. H. Kinetics of decomposition of Colorado oil shale. II. Carbonate minerals. LLNL Report UCRL-52089 Part 2; Lawrence Livermore National Laboratory (LLNL): Livermore, CA, 1978.
253
Energy Fuels 2011, 25, 251–259
: DOI:10.1021/ef101225g
Bauman and Deo
Table 1. High and Low Molecular Weights and Associated Stoichiometry for Reactions 1 and 2
molecular weight (þ) molecular weight (-) formula (þ) formula (-) stoichiometric reaction 1 (þ) stoichiometric reaction 1 (-) stoichiometric reaction 2 (þ) stoichiometric reaction 2 (-)
kerogen
heavy oil
light oil
gas
20000.55 2974.84 C1479H2220 C220H330 -1 -1
424.49 424.61 C31.75H42.82 C31.76H42.81 37.29 5.55 -1 -1
152.03 151.99 C11.19H17.51 C11.19H17.50 13.86 2.06 2.18 2.18
52.01 51.95 C3.35H11.63 C3.35H11.62 25.03 3.72 0.06 0.06
been used in a wide variety of applications.20,21 These experimental design tools have also been applied to various oil reservoir studies.22 A common experimental design is the 2k full factorial design. These designs test k factors at two levels for each factor: high and low levels. Each combination of high and low values of each factor is called a run. Full factorial designs require 2k runs to test every possible combination of high and low levels for k factors. When the number of factors is excessive or runs are expensive, fractional factorial designs are used for efficiency. In fractional factorial designs, runs are selectively eliminated with the assumption that higher order interactions are much less significant than individual factors without interactions. These designs are represented as 2k-p fractional factorial designs. Fewer runs are required, but information about the significance of higher order interactions is confounded with information about individual parameters. Experimental design and analysis methods are useful for comparing the sensitivity of a response because of variable input parameters, including their possibly significant interactions. The parameters of particular interest in this study are molecular weight of kerogen, activation energy for kerogen cracking reaction 1, activation energy for heavy oil cracking reaction 2, activation energy for light oil cracking reaction 3, activation energy for gas cracking reaction 4, distributed representation for activation energy for kerogen cracking reaction 1, relative permeability representation, and reaction enthalpy. Each of these parameters is required for calculating the mass, energy, and momentum balances solved by the simulator. Activation energies are required for calculating the reaction rate term in the mass balance equation. Relative permeability is used to calculate the flow term with Darcy’s law. Reaction enthalpy is incorporated in the energy balance equation. Ranges for each of these parameters were estimated from various literature data, inherent uncertainty, or are estimated to explore sensitivities. The molecular structure of kerogen is largely unknown. The molecular weight of kerogen has been reported in ranges from about 300023 to 27 000.24 The stoichiometry in the chemical reactions and the initial concentration of kerogen are dependent upon the choice for molecular weight of kerogen to
methane 16.04 16.04 CH4 CH4 17.06 2.54 0.03 0.03
char
coke
12.60 12.55 CH0.55 CH0.53 38.71 5.8 7.13 7.13
14.55 14.55 C1.19H0.32 C1.19H0.32 0 0 0 0
Table 2. Range of Activation Energies reaction
low activation energy (kJ/mol)
high activation energy (kJ/mol)
kerogen cracking (1) heavy oil cracking (2) light oil cracking (3) gas cracking (4)
195 208 208 235
225 260 260/233 270
conserve volume and mass for all simulation runs. Consequently, when the molecular weight of kerogen is changed between simulation runs, the stoichiometry of the reactions and the initial molar concentration of kerogen in the pore space must also be adjusted for mass and volume consistency. The range of molecular weight and associated stoichiometry for reactions 1 and 2 are shown in Table 1. The values in Table 1 demonstrate that stoichiometry for such a kinetic mechanism depends upon the properties (molecular weight and H/C ratio) of the pseudo-components and is therefore “non-unique”. Mass and elemental balances are important in these representations. Ranges for appropriate activation energies have been reported18 and can vary significantly depending upon experimental methods and analysis of results. These ranges for activation energy are shown in Table 2. Studies have reported that activation energy for kerogen pyrolysis is most appropriately modeled with some distribution,18,25 but it is uncertain how much impact different representations of activation energy have on the simulation results at large scales. A normal distribution with 5 kJ/mol standard deviation is shown in Figure 2. Distribution of activation energies for kerogen pyrolysis is a complex function but is sometimes represented by the normal distribution.26 A perfect activation energy distribution cannot be represented exactly in STARS; therefore, discrete quantities, determined by integrating under the distribution curve, represent kerogen reacting with a specified activation energy according to the distribution. Relative permeability representations are often approximated in simulation, but such approximations may have significant implications. The range of relative permeability curves are shown in Figure 3, with the low level being more linear and the high level being curved. The shape of the relative permeability curves depends upon the resource, the wetting characteristics of the rock, and the constituents present in the pore space. Finally, the heat of reaction could play an important role in the heat-transfer efficiency depending upon the characteristics of the associated reactions. Efficient heat transfer through an oil shale reservoir is crucial to any successful operation. Heat of reaction for oil shale pyrolysis
(20) Box, G. E. P.; Hunter, W. G.; Hunter, J. S. Statistics for Experimenters: An Introduction to Design, Data Analysis, and Model Building; John Wiley and Sons, Inc., New York, 1978. (21) Lawson, J.; Erjavec, J. Modern Statistics for Engineering and Quality Improvement; Duxbury Press: Pacific Grove, CA, 2001; pp 465-513. (22) Peng, C. Y.; Gupta, R. Experimental design and analysis methods in multiple deterministic modelling for quantifying hydrocarbon inplace probability distribution curve. Proceedings of the Society of Petroleum Engineers (SPE) Asia Pacific Conference; Kuala Lumpur, Malaysia, March 29-30, 2004; SPE Paper 87002. (23) Yen, T. F.; Chilingar, G. V. Introduction to oil shales. In Oil Shale; Yen, T. F., Chilingar, G. V., Eds; Elsevier Science Publishing Company: Amsterdam, The Netherlands, 1976; pp 181-198. (24) Behar, F.; Vandenbroucke, M. Chemical modelling of kerogens. Org. Geochem. 1987, 11, 15.
(25) Sundararaman, P.; Merz, P. H.; Mann, R. G. Determination of kerogen activation energy distribution. Energy Fuels 1992, 6 (6), 793– 803. (26) Burnham, A. K.; Braun, R. L. Global kinetic analysis of complex materials. Energy Fuels 1999, 13 (1), 1–22.
254
Energy Fuels 2011, 25, 251–259
: DOI:10.1021/ef101225g
Bauman and Deo
Figure 2. Normal distribution of activation energies for reaction 1.
Figure 3. Oil/water and liquid/gas relative permeability curves. Sw is water saturation, and Sl is liquid saturation. Gas (krg), water (krw), oil/ water (krow), and oil/gas (krog) relative permeability parameters are used in Stone’s II relative permeability model.
has been reported,27 but it is not certain how much heat is “lost” to the reaction compared to heat required to raise the rock temperature. Results and Discussion
row represents a simulation run at the parameter levels specified in the run. The uncoded levels for these parameters are shown in Table 3. Details about each of these parameters have been described. The numerical performance for each of these runs was not equal. Some of the runs had
The initial experimental design was a 27-4 fractional factorial design. The eight run design for the initial 7 factors [excluding (8) heat of reaction] is shown in Figure 4. Each
(27) Camp, D. W. Oil shale heat capacity relations and heats of pyrolysis and dehydration. Proceedings of the 20th Oil Shale Symposium; Golden, CO, 1987.
255
Energy Fuels 2011, 25, 251–259
: DOI:10.1021/ef101225g
Bauman and Deo
Figure 4. Eight run fractional factorial screening experimental design. Table 3. Uncoded Parameters for Screening Design factor
physical parameter
X1
molecular weight/ stoichiometry/concentration Eact reaction 1 Eact reaction 2 Eact reaction 3 Eact reaction 4 Eact distribution reaction 1 relative permeability
X2 X3 X4 X5 X6 X7
low level (-)
high level (þ)
3000
20000
195 208 208 235 without linear
225 260 260 270 with curved
Figure 6. Aerial view of the simulated wedge. The distance between heating wells was reduced to 26.5 ft.
Figure 7. Normal probability plot of the effects on the ultimate recovery of oil from 27-4 fractional factorial design.
Figure 5. Pareto chart for parameter effects on simulation time.
excessive time-step cuts because of rapid changes in gas saturation. The response chosen for these runs was simulation time to pinpoint the possible causes of these time-step reductions. Figure 5 is a Pareto chart displaying the impact that each of these parameters has on the simulation time. Activation energy for reaction 3 or factor X4 had the greatest impact on the simulation time. After investigation, it appeared that the simulation time increased significantly when the activation energy for reaction 3 was greater than the activation energy for reaction 4. This could be due to the combination of rapid gas creation coupled with high gas mobility, causing rapid gas saturation changes. The high value for factor X4 was lowered to 233 kJ/mol, as shown in Table 2, and no major differences in simulation time were observed in subsequent runs. Using the lowered value for factor X4, runs in the fractional factorial design were completed with ultimate recovery of oil as the output response. None of the simulations produced acceptable amounts of oil. Upon inspection, it was found that oil generated from kerogen had inadequate mobility in lower temperature zones far from the heaters to flow to the producer. As a result, oil components had large residence times in the reservoir and eventually converted further to gas and residual components. This result gives insight into the design of such a process, specifically the spacing needed between wells for successful operation. If heating wells are drilled too
far from producing wells, the residence time of the oil in hot zones of the reservoir will be excessive and these oils will convert to gases or residual solids, significantly reducing or even prohibiting the production of oil. However, capital and operating costs increase with the number of wells drilled. Well spacing is a crucial design consideration for this process because excessive residence time of products in the reservoir and the cost of drilling wells are competing considerations for optimal process design. The initial dimensions of the simulated domain were changed to resolve this issue of excessive residence time of the oil. Reducing the spacing between these wells assures that the whole reservoir was at a high enough temperature for adequate oil mobility. Figure 6 shows the modified dimension of the simulated wedge, changing the distance between heaters from 53 to 26.5 ft. The same fractional factorial design was used with ultimate recovery of oil as the response. The normal probability plot in Figure 7 illustrates the results of the runs. Normal probability plots, like Pareto charts, are useful by visualizing the significance of the effects for each factor. Dominating effects will appear as outlier points on a normal probability plot. The points of the effects on this plot in Figure 7 are linear without outliers, indicating that there is no evidence from these runs that any factors are dominant or insignificant. With the 27-4 fractional factorial design used, single-factor effects are confounded with pair interaction 256
Energy Fuels 2011, 25, 251–259
: DOI:10.1021/ef101225g
Bauman and Deo
Figure 8. Experimental design for 6-8 factors without confounding of individual parameters. Table 4. Summary of Calculated Effects
Figure 9. Pareto chart from 16 run fractional factorial design for 8 factors.
effects and higher order interactions. Additional runs are necessary to isolate the effects of individual parameters from confounding with the effects of higher order interactions. Further runs were performed with a 16 run fractional factorial design for 6-8 factors. All 8 factors, including heat of reaction, were tested with this design. The design used is shown in Figure 8, where factors E1-E7 represent possible interactions between parameters but individual factors are isolated from possible confounding with interactions. For example, E1 in Figure 8 could be significant because of any combination of interactions between X1 and X2, X3 and X5, X4 and X8, or X6 and X7. The results from these runs are displayed in a Pareto chart in Figure 9. It appears that the most significant factors are X2, X4, X6, and X7 (please see Table 3 for uncoded parameters), along with higher order interactions between parameters, likely between these most significant factors. These factors are activation energy for reaction 1, activation energy for reaction 3, activation energy distribution representation for kerogen conversion, and relative permeability representation. It appears that activation energy for reaction 4 and heat of reaction have the least impact on ultimate recovery of oil. Reaction 4 is important in determining product distributions, but for the dynamics of the reaction and flow considered in these simulations, the oil recovery (main outcome) is not affected by the variability in the activation energy of reaction 4. The results from these runs can be used in a 24 full factorial design without any additional runs. The data were regressed
factors
β
intercept X1 X2 X3 X4 X1X2 X1X3 X1X4 X2X3 X2X4 X3X4 X1X2X3 X1X2X4 X1X3X4 X2X3X4 X1X2X3X4
294.6071 -25.9846 48.71381 -8.88369 -36.7407 17.66644 -19.9203 7.082688 14.89081 -2.82519 1.098062 -5.35956 0.575937 -3.37056 -0.70869 2.191438
with the polynomial model shown in the equation below, where β0 is the intercept (global mean), β is single and higher order interaction linear coefficients, and x is input variables. This polynomial model forms a multivariate surface called a response surface. The effects are calculated by taking the difference of the averages of the responses at high and low levels of each factor, and at high and low levels for interacting factors. The coefficients β are half of the calculated effects. The coefficients are summarized in Table 4. k k X X X βi, j xi xj βi xi þ y ¼ β0 þ i¼1
þ
XX X
i6¼j6¼l
i6¼j
βi, j, l xi xj xl þ βi, j, l, m xi xj xl xm
This model fit the experimental output data exactly. Although this is not a theoretical model and may have little physical significance, insight about the significance of each parameter in the explored ranges can be garnered. Typically, higher order linear interaction effects are assumed to be negligible and can be used to estimate error.22 Expert opinion and knowledge are advantageous for estimating error, and elimination of terms in this model perhaps is not justified because this knowledge is unknown.22 Three random validation simulations within the experimental space were run to estimate the quality of the response surface, the empirical regression model, compared to a 257
Energy Fuels 2011, 25, 251–259
: DOI:10.1021/ef101225g
Bauman and Deo
Figure 10. Histogram of Monte Carlo calculations of the response surface.
STARS simulation. The difference between the response surface approximations for ultimate oil recovery and STARS simulation results ranged from 3 to 15%. The quality of the response surface could be improved at the cost of more experimental runs, by either reducing the experimental space or adding additional runs to estimate curvature because of nonlinearities when parameters are continuous. Alternative experimental designs could possibly provide more accurate response surfaces with comparable or fewer total runs; however, many of these alternative designs require additional expert knowledge about the problem or unjustified assumptions. Monte Carlo simulations were performed to characterize the response surface. Random values for each parameter with uniform distributions were chosen for each run. A histogram of 80 000 Monte Carlo runs is shown in Figure 10. The average value in these runs was 294.7 barrels of oil, with a standard deviation of 40.1 barrels of oil. A normal distribution with these values is also shown in the figure for comparison. It appears that the Monte Carlo results are slightly skewed to the right of a normal distribution. This exercise helps to quantify the effects of variations in input parameters on the desired output. The shape of this distribution could be affected by the response surface itself, the sampling locations for Monte Carlo simulation, or the distributions assigned to each of the factors.
distribution representation have significant impacts on the ultimate recovery of oil in these simulations. Expert knowledge or similar studies including large-scale physical experiments are important for estimating statistical error for developing validated surrogate models. Otherwise, more runs are necessary for improving the quality of these models for approximating simulator results. The interplay between various flow and kinetic parameters has been explored. Geomechanical, heat-transfer, and equilibrium parameters, for example, may also play significant roles at certain scales in production results for such complex reactive transport systems. Parameters from acceptable theoretical models can also be included in experimental designs to evaluate their impact on results and to include these parameters in constructing response surface approximations, as illustrated in this paper. Response surfaces can be characterized to quantify risk and uncertainty of simulations according to variation in input data. Acknowledgment. The authors acknowledge financial support from the National Energy Technology Laboratory, U.S. Department of Energy (Grant DE-FE0001243). The authors also thank Computer Modeling Group Limited (CMGL), Calgary, Alberta, Canada, for providing academic licenses to their reservoir simulators.
Nomenclature Conclusions
Eact = activation energy γp = phase-specific gravity ^ p = phase enthalpy H ˆ ΔHr = heat of reaction K = permeability tensor krp = phase relative permeability μp = phase viscosity Nf = number of fluid phases Np = number of phases Nr = number of reactions Pp = phase pressure p = phase index φ = porosity Qe = source term for energy Qm,i = source term for mass of component i
Although results for oil shale simulations in this study are calculated with theoretical governing equations, the interplay within various parameters is not trivial because of competing physical phenomena. Combinations of parameters that expose possible competing phenomena can have significant numerical implications. Molecular representations for kerogen with associated stoichiometry, heat of reaction for kerogen decomposition, intermediate oil cracking reaction (reaction 2) activation energy, and continuing gas cracking (reaction 4) reaction activation energy are insignificant in determining the ultimate recovery of oil at the scale simulated in this paper. Kerogen cracking (reaction 1) activation energy, relative permeability representation, oil cracking to gas (reaction 3) activation energy, and activation energy 258
Energy Fuels 2011, 25, 251–259
: DOI:10.1021/ef101225g
Bauman and Deo
T = temperature t = time Up = phase internal energy Urock = rock internal energy νBp = phase velocity xp,i = mole fraction of component i in phase p z = height or depth
RCi = residual mass for component i RE = residual energy R = reaction rate r = reaction index Fp = phase molar density Sp = phase saturation Sr,i = stoichiometric factor of component i in reaction r
259
