Article pubs.acs.org/EF
Hybrid Multiscale Modeling through Direct Substitution of PoreScale Models into Near-Well Reservoir Simulators Tie Sun, Yashar Mehmani, and Matthew T. Balhoff* Department of Petroleum and Geosystems Engineering, The University of Texas at Austin, 1 University Station, C0300, Austin, Texas 78712-0228, United States ABSTRACT: A hybrid, multiscale reservoir simulator has been developed that directly substitutes pore-scale models for continuum-scale, finite-difference grids. Over 7500 pore scale models (∼75 million pores) are inserted into a ∼1 m2 pore-scale region around a producing well, which is coupled to an outer, fine-scale finite-difference model with 10 000 cell-centered grids. This computationally demanding problem is possible because the domain is decomposed into manageable subdomains (thousands of pore-scale models and one continuum region) and then coupled together using finite-element mortars to ensure continuity between the subdomains. Here, flow is single-phase, single-component, and steady-state, but the potential of the simulator for more complex, near-well flows is demonstrated. A priori, direct upscaling and a posteriori, global upscaling are both performed, and it is shown that the direct approach results in significant inaccuracies in the pressure and flow fields but the global upscaling approach results in a good approximation.
1. INTRODUCTION Sub-surface simulators are used to model flow and transport of many applications in porous media, including hydrocarbon recovery, carbon sequestration, groundwater remediation, and nuclear waste storage. These simulators are modeled at the continuum scale using semi-empirical equations (e.g., Darcy’s law) because of the large scales (103−104 m) of interest in the sub-surface. The models require macroscopic inputs, such as permeability, porosity, and relative permeability curves. In general, these values must be estimated from the geological data or measured from laboratory experiments. Recent improvements in pore-scale (10−4−10−2 m) modeling suggest that direct numerical simulation through the pore space has reached a level of accuracy that they could be used as surrogates for experiments.1−5 This is largely a result of improved imaging capabilities and better description and solution of the fundamental flow equations in the pore space. Pore-scale approaches include direct computational fluid dynamics modeling using the finite-element method (FEM) or Lattice Boltzmann,6 level-set approaches,7 and network modeling.8,9 Network models require some approximations to the flow equations but have the advantage of modeling at slightly larger scales (10−3−10−2 m). Despite improvements in pore-scale modeling, questions remain as to the ability to extract macroscopic parameters from a stand-alone, small-scale model. Even if the model is large enough to be considered a representative elementary volume (REV), inaccuracies can arise because the imposed boundary conditions do not reflect flow and transport in surrounding media.10−13 Ideally, flow and transport could be modeled in a single, large pore-scale model that spans 10−2−100 m and millions of pores, but this is computationally impractical. As a result, accurate multiscale methods are sought that upscale pore-level flow behavior to the continuum scale. A lofty goal would involve multiscale simulation via direct substitution of pore-scale models for fine grids in strategic parts of a reservoir (e.g., near wells). Such a hybrid (pore-to-continuum) model © 2012 American Chemical Society
would be unique and allow for extremely powerful modeling of sub-surface flow phenomena. The difficulty of such an approach includes computational limitations but, more importantly, an ability to bridge and couple scales. Balhoff et al.10 were one of the first to couple a pore-scale model directly to a continuum model. That work demonstrated that boundary conditions provided by adjacent media do affect pore-scale flow. Scheibe et al.14 coined the term “hybrid methods” to refer to coupled modeling of both scales in the context of porous media simulation; a research area still in its infancy. Tartakovsky et al. 15 used a smoothed particle hydrodynamics (SPH) formulation of diffusion-dominated reactive transport in porenetwork and continuum models and managed to couple them explicitly without iterations. For a review of particle-based multiscale simulations, see the study by Koumoutsakos.16 Chu et al.17 used a novel way of coupling mass conservation equations with constitutive laws derived directly from porescale networks. Their work is based on sampling small-scale representations of the domain. Proofs of convergence were provided, and examples for case studies were given. Hybrid multiscale modeling can also be accomplished using domain decomposition approaches, as was performed by Balhoff et al.10 Mortar coupling is one class of domain decomposition methods in which a large boundary value problem (BVP) is split into a number of smaller BVPs defined over smaller non-overlapping subdomains, with interface conditions in between them. The entire problem along with the interface conditions is solved through an iterative procedure in an effort to match fluxes at the interfaces. The decoupled subdomains are independent from each other; therefore, they can include different models, physics, solution algorithms, or meshes. Moreover, the subdomains can have non-matching Received: June 13, 2012 Revised: July 19, 2012 Published: July 27, 2012 5828
dx.doi.org/10.1021/ef301003b | Energy Fuels 2012, 26, 5828−5836
Energy & Fuels
Article
dimensions. This allows for immense flexibility in modeling highly heterogeneous multiprocess phenomena in the subsurface with the ability to increase modeling resolution locally. Most importantly, the problem easily lends itself to parallel computing because of the independence of the subdomains. Arbogast et al.18 used mortar spaces to approximate the trace of the solution at the subdomain interfaces by imposing weak flux continuity. In this study, they set up a simple elliptic model problem defined on non-matching multiblock grids and used it to derive the weak formulation of the subdomains as well as the interface problem and investigated convergence and error estimates. Wheeler et al.19 and Peszynska et al.20 demonstrated the applicability and capability of the mortar method in various multiphase case studies in highly heterogeneous domains, demonstrated its potential as a new generation reservoir simulator (IPARS), and denoted them as a new upscaling methodology for multiphase flow, with similarities to the “subgrid upscaling” method. 21 Thomas and Wheeler 22 introduced a novel extension of an enhanced velocity mixed FEM to flow and reactive transport problems. Balhoff et al.11 first used mortar coupling to couple pore-scale models with other pore-scale models or continuum models. In that work, pore-scale models were coupled to adjacent finitedifference models and interface pressure fields were determined in such a way that fluxes matched weakly over that same interface. Sun et al.13 used the mortar approach to couple hundreds of pore-scale models in a rigorous upscaling scheme. They showed that the upscaled permeability was more accurate than obtained through more straightforward hierarchical methods. However, to our knowledge, no one has successfully attempted to directly substitute pore-scale models into a subsurface reservoir simulator, at least at the scale studied here. In this work, we present a first-of-its-kind sub-surface, hybrid reservoir simulator that directly substitutes thousands of porescale models for continuum-scale grids in a ∼1 m2 region. The pore-scale models are coupled to each other as well as to an outer Darcy-scale region using mortars. The novel multiscale simulator could be used for accurate modeling of nonlinear and dynamic near-well phenomena (non-Darcy flow, skin, acid transport in injection strategies, formation damage, etc.) that require pore-scale investigation while still maintaining the computational efficiency of continuum-scale modeling in the remainder of the reservoir.
Figure 1. Areal geometry and boundary of the hybrid, multiscale simulator. The continuum region is 2D with a thickness of 1 cm, and the pore-scale region is 3D. continuum-scale region includes 10 000 fine (1 cm3) grid blocks of uniform permeability (1 × 10−6 cm2). Constant pressure boundary conditions are imposed in the well (P = 7 × 105 Pa at r = 10 cm) and the outer boundary (P = 7.2 × 105 Pa). The simulator is decomposed into 7549 subdomains (7548 porescale network models and 1 continuum-scale, finite-difference model). Each of the subdomains is solved independently and could be distributed across multiple processors for computational efficiency, although here, they are solved on a single processor. Flow is singlephase and steady-state throughout the model. Each network model is solved in the usual way;1,8 a mass balance is imposed at each of the N pores in the network, which leads to a system of N linear equations. Solution to the system of equations gives pore pressures, which, in turn, provide throat flow rates and total Darcy velocity across all six boundaries of the network. The continuum, finite-difference region is solved implicitly using cell-centered grids, Darcy’s law for velocity, and a harmonic mean for interblock permeability. Dirichlet boundary conditions are imposed at four subdomain interfaces, and Neumann boundary conditions are applied at the top and bottom sides. Unlike typical network modeling studies, the boundary conditions at network faces are not a single, constant pressure but, instead, are pressure fields that vary spatially over the network face. The pressure boundary conditions are described as a linear combination of basis functions over a coarse mortar grid. Here, we have used 3 × 3 mortars at network interfaces and implemented piecewise, bilinear basis functions
2. MATHEMATICAL APPROACH 2.1. Multiscale Simulator Description and Modeling Approach. Figure 1 depicts the geometry of the near-well, hybrid, and multiscale simulator, which is areal and pseudo-three-dimensional (3D) [the pore-scale models are 3D, but the continuum blocks are only 1 grid block in thickness and therefore two-dimensional (2D)]. Although the geometry is Cartesian and square, the wellbore geometry is a multisided polygon, so that it approximates a circle with a radius of 10 cm. The pore-scale region extends to a radius of 0.5 m, and the continuum region extends an additional 0.5 m. The pore-scale region includes 7548 pore-scale, network models, each with dimensions of 1 × 1 × 1 cm. Therefore, the pore-scale region covers an area of 0.75 m2. Although 7548 unique network models could have been used, for simplicity, 40 different networks were repeated randomly in the domain. All but one of the networks were mapped from computergenerated sphere packings using a Delaunay tessellation,23 and the last network was mapped from a real sandstone,24 imaged using X-ray microtomography. This sandstone network added heterogeneity to the problem. The network models ranged from 400 to 40 000 pores, and permeability ranged from 1 × 10−7 to 9 × 10−5 cm2. In all, the porescale region in the simulations contains 30−75 million pores. The
Ndof
pi =
∑ αjϕj j=1
(1)
where p is the spatially varying pressure, Ndof is the number of degrees of freedom, and ϕ is the piecewise basis function. Defining the constants (α values) gives the pressure field over the entire network face. The pressure field is then projected onto the network model; pressure at the boundary pores at given spatial locations can be calculated from the boundary pressure field (eq 1). The network model can then be solved using these boundary pore pressures. The multiscale simulator is able to solve subdomains independently and efficiently while ensuring continuity of pressures and fluxes between the subdomains. Mortars are coarse finite-element spaces that glue the subdomains together to ensure continuity, in a weak finiteelement sense. The mathematical details are summarized in several references,18,20,11 but to summarize, the coefficients of the basis functions in eq 1 are guessed, the subdomains are solved using this 5829
dx.doi.org/10.1021/ef301003b | Energy Fuels 2012, 26, 5828−5836
Energy & Fuels
Article
boundary condition, the jump in fluxes across the mortar space is computed, and then the coefficients are iterated upon until the flux jump is 0 (or below a predetermined tolerance). For linear problems, such as this one, only one iteration is required. Mortars have been successfully used to couple continuum models to continuum models,18,20 pore-scale models to pore-scale models,11,13 and porescale models to continuum models.11 The advantage of the mortars here is that one can model a pore-scale region orders-of-magnitude larger (30−75 million pores and 0.75 m2) than previously investigated and directly substitute and couple it with a conventional, finitedifference reservoir simulator. 2.2. Upscaling. The multiscale simulator is unique and has the capability to capture complex heterogeneity and model unique flow and transport physics at the pore scale without the loss of information via upscaling. Nonetheless, the simulation is computationally demanding (the pore-scale region alone has up to 75 million unknowns). Ultimately, upscaling may need to be performed, so that simulations can be conducted more efficiently. We propose two types of upscaling: a priori, direct upscaling and a posteriori, global upscaling. The two approaches are compared to the full, multiscale solution. Upscaling of pore-scale models is usually performed in the following, direct fashion: a constant pressure gradient is imposed in one direction; no-flow or periodic boundary conditions are imposed in the other two directions; and macroscopic properties (e.g., permeability) are back-calculated from the resulting Darcy velocity using a macroscopic model (e.g., Darcy’s law). Here, we have performed these calculations for all 40 networks in two directions; therefore, directionally dependent, diagonal permeability tensors are obtained and then inserted into the finite-difference reservoir grid in place of the network models. Note that flow effects of surrounding media are not included in the boundary conditions because no coupling between models is used. This a priori, direct upscaling is the usual approach to extract macroscopic properties from pore-scale simulations.2,4,9,25 A second upscaling approach is employed and investigated here, one that involves a posteriori upscaling and better accounts for the impact of surrounding media. The a posteriori upscaling is a global upscaling approach.26 That is, once the full, multiscale simulator is solved, the permeability of each individual network model is backcalculated from the resulting flow and pressure boundary conditions
Kx , ij =
K y , ij =
solution for steady-state, radial flow around a well. The wellknown solution is given as P(r ) = Pw + (P0 − Pw )
Figure 2. Comparison between hybrid, multiscale solution using uniform, homogeneous, and isotropic networks and the analytical solution for radial flow (eq 4).
(2)
3.2. Homogeneous Pore-Scale Region. Simulations (case II) were performed for a relatively homogeneous porescale region, but with different networks scattered throughout. All of the networks used were mapped from computer-based, unconsolidated sphere packings, and the sandstone was not used. The models may represent loose sand or granular media. The networks still exhibit some pore-scale heterogeneity but are significantly more homogeneous than the sandstone network. The pressure field obtained from the multiscale, near-well model from simulation case II is presented in Figure 3. The colored points in the near-well region represent pore pressures from the pore-scale models (size of the points is arbitrary, and white space is shown because pores are discrete), and the pressures in the outer region represent pressure values of continuum grids. The pressure and flux fields obtained from the multiscale model are compared to the two upscaling approaches (direct upscaling of the finite-difference model and the a posteriori, global upscaling model). The comparison of the pressure (Figure 3) and flow fields (Figure 4) suggests that the direct upscaling approach produces a good approximation to the multiscale simulation. This is further verified in the quantitative comparison of various
qy , ijμΔy AΔpy , ij
(4)
where P(r) is the pressure at some radius r, Pw and P0 are the boundary pressures at the well and outer boundary, respectively, and rw and r0 represent the radius of the well and outer boundary, respectively. The hybrid solution was performed using 7548 identical, relatively homogeneous, and isotropic pore-networks (case I) derived from uniform sphere packings (direct upscaled permeability = 1.05 × 10−5 cm2). The continuum-scale region had the same permeability. As shown in Figure 2, the hybrid solution matches very well with the analytical solution. The relative error (L2 norm) between the analytical and hybrid pressure solution is 0.04%.
qx , ijμΔx AΔpx , ij
log(r /rw ) log(r0/rw )
(3)
where qx,ij and Δpx,ij are the flow and pressure drop along the x direction extracted from the multiscale model for the ijth network model. The boundary pressures (used to compute Δp) on the network were computed by taking a weighted average of the nodal pressures on the mortar interface grid. In contrast to direct upscaling, where Kx and Ky are estimated a priori for the 40 repeated network models, 7548 sets of Kx and Ky are estimated from the a posteriori, global upscaling approach.
3. RESULTS AND DISCUSSION The pressure and flux fields obtained from the multiscale simulator are compared to the two upscaling approaches (a priori, direct upscaling and a posteriori, global upscaling). The results are presented in the form of pressure contours, quantitative pressure profiles along four major angles, flow contours, and total production rates (flow rate into the well). 3.1. Validation with Analytical Solution. The hybrid, multiscale simulation was first verified against the analytical 5830
dx.doi.org/10.1021/ef301003b | Energy Fuels 2012, 26, 5828−5836
Energy & Fuels
Article
Figure 3. Pressure profile of simulation case II derived from three approaches: (a) multiscale, (b) direct upscaling, and (c) a posteriori, global upscaling. Note that, in panel a, the points represent pore pressure in the pore-scale models. The pressure unit is Pa.
Figure 4. Flow field of simulation case II derived from three approaches: (a) multiscale, (b) direct upscaling, and (c) a posteriori, global upscaling. Units are m3/s.
Figure 5. Pressure profile for simulation case II at four angles around the well, 0°, 90°, 180°, and 270°, from the positive x axis extracted from the three models. The solid line indicates the direct upscaling model; the dashed line indicates the a posteriori, global upscaling model; and the black points represent the multiscale model. Note that, in this simulation, pressures from all three models closely match each other. The pressure unit is Pa.
5831
dx.doi.org/10.1021/ef301003b | Energy Fuels 2012, 26, 5828−5836
Energy & Fuels
Article
Figure 6. Pressure profile of simulation case III derived from three approaches: (a) multiscale, (b) direct upscaling, and (c) a posteriori, global upscaling. Note that the direct upscaling model yields more pressure drop in the pore-scale subdomain compared to the multiscale model and global upscaling model. The pressure unit is Pa.
Figure 7. Flow field of simulation case III derived from three approaches: (a) multiscale, (b) direct upscaling, and (c) a posteriori, global upscaling. Units are m3/s.
Figure 8. Flow rate difference between (a) the multiscale simulator and the direct upscaling approach and (b) the multiscale simulator and the a posteriori, global upscaling approach in case III. Note that the scale is logarithmic. Units are m3/s.
where stronger pore-scale heterogeneity is incorporated, the direct upscaling result shows strong discrepancy with respect to the multiscale simulation. 3.3. Heterogeneous Pore-Scale Region. The last simulations (case III) were performed by using the much more heterogeneous sandstone. It was randomly repeated 4500 times (60% of the networks), and the remaining 39 sphere-pack networks comprised the remaining 3000 pore-scale subdomains. Figures 6 and 7 show the pressure and flux, respectively, around the well obtained from the three simulation approaches.
pressure profiles (Figure 5) at four angles from the horizontal (0°, 90°, 270°, and 360°). This observation suggests that direct upscaling can be accurate for single-phase, Newtonian flow if the pore-scale models are relatively homogeneous and isotropic. The a posteriori, global upscaling approach also successfully reproduces the multiscale model results. The advantage of upscaling is that, once upscaled permeabilities are computed, continuum simulations can be performed, which are much less computationally expensive than the multiscale model. Quantitatively, both upscaling methods resulted in less than 1% error in well production rates. However, in later examples, 5832
dx.doi.org/10.1021/ef301003b | Energy Fuels 2012, 26, 5828−5836
Energy & Fuels
Article
Figure 9. Pressure profile for simulation case III at four angles around the well, 0°, 90°, 180°, and 270°, from the positive x axis extracted from the three models. The solid line indicates the direct upscaling model; the dashed line indicates the a posteriori, global upscaling model; and the black points represent the multiscale model. The pressure unit is Pa.
Figure 10. Permeabilities (cm2) of each network obtained from the two upscaling approaches for case III. The red triangles represent direct upscaling permeability, and the blue dots represent permeabilities obtained from the a posteriori, global upscaling approach.
permeability for the pore-scale region. The flow paths, however, are captured well by the global upscaling approach (Figure 7c), and nearly identical flux fields are depicted. This is further confirmed by Figure 8, which depicts the error (difference in flux for multiscale solution and upscaled solution). Quantitatively, the a posteriori method underestimates production by only 4%, while the a priori method did so by over 40%. The pressure profile (Figure 9) obtained by the a posteriori, global upscaling approach (blue lines) closely matches at the four selected angles (0°, 90°, 270°, and 360°) to the multiscale simulation (black points), suggesting that global upscaling accurately reproduces the multiscale model results. The upscaled permeabilities for case III are summarized in Figure 10. In this figure, the x axis represents the network index
The a posteriori, global upscaling approach (Figure 6c) produces similar results as the multiscale model (Figures 6a). However, the direct upscaling approach (Figure 6b) overestimated the pressure drop in the pore-scale region (permeability, on average, is underestimated). The simplified boundary conditions employed may be significantly different from those obtained by coupling to surrounding media. Figure 7 provides the flow contour plots obtained from the three approaches. A few high flow paths are observed in the pore-scale region for the full, multiscale simulation (Figure 7a). They can also be observed in the direct upscaled solution (Figure 7b), but the flow values are significantly smaller than those in the multiscale simulation. This observation also confirms that the direct upscaling underestimated the 5833
dx.doi.org/10.1021/ef301003b | Energy Fuels 2012, 26, 5828−5836
Energy & Fuels
Article
Figure 11. Histograms for 2 of the 40 networks used in the simulations (one homogeneous, sphere pack network and the other, the sandstone network) for both the x and y directions. Panels a and b show the permeability histogram for the heterogeneous sandstone network. Panels c and d show the permeability histogram for a relative homogeneous uniform sphere pack network. The histograms are overlaid with direct upscaled permeability (indicated by the dark line) and the geometric mean of the globally upscaled permeability (indicated by the gray line).
well region, preferential flow paths will form as wormholes are created in the pore space. This is a pore-scale phenomenon by nature; proper modeling requires reaction with the grain space to form expanded pore channels.27 The prediction of wormhole growth and resulting permeability improvement would be difficult (or impossible) to model at the continuum scale. NonDarcy, Forchheimer flow is another near-well application that may require detailed pore-scale modeling. Some studies4 suggest that a macroscopic model may not suffice for modeling these high-velocity flows. Injection of CO2 in sequestration strategies and hydrocarbon production in heterogeneous, damaged formations may also benefit from a hybrid simulator. Mortars are a critical component to the feasibility of this hybrid simulator. It is unlikely that a pore-network model could be created, much less solved, for the size investigated here (tens of millions of pores and nearly a square meter in area). The mortar method allows for the pore-scale region to be decomposed and solved in small, manageable parts. This is an advantage even on a single processor but would be most efficient in a parallel computing environment because this is nearly an embarrassingly parallel problem. The coarse, interface problem used to ensure continuity is fast compared to the subdomain problems; therefore, each subdomain (or a group of subdomains) could be sent to an individual processor. The simulator is not without limitations. Here, the pore-scale region is 0.75 m2 in area, but the pore-scale models are
and the y axis represents the effective, upscaled permeability values. Because the a posteriori approach incorporates the impact of surrounding media, multiple permeability values are obtained for the same network at different locations in the reservoir. Rather than one single value, the globally-upscaled permeabilities span over 2 orders of magnitude for most of the network models, which also verifies that an accurate upscaled permeability may be significantly affected by the neighbors/ boundary conditions in addition to the pore-scale model itself. In particular, when heterogeneous network samples are present at complex flow conditions (i.e., nonlinear dominant flow region), the permeabilities vary significantly. Figure 11 shows histograms for 2 of the 40 networks used in the simulations (one homogeneous, sphere pack network and the other, the sandstone network) for both the x and y directions. The direct, upscaled permeability matches well with the geometric mean of the global, upscaled permeability for the more homogeneous network. However, a larger discrepancy is observed for the heterogeneous (sandstone) network sample, which partially explains the error obtained using direct upscaling. 3.4. Discussion of Applications. The hybrid, multiscale simulator developed here is demonstrated for simple physics (single-phase, single-component flow). The real use of this type of simulator will likely be for more nonlinear and dynamic flows. For example, if acids are injected to stimulate the near5834
dx.doi.org/10.1021/ef301003b | Energy Fuels 2012, 26, 5828−5836
Energy & Fuels
■
relatively large (1 cm3 as opposed to 1 mm3). Likewise, the continuum blocks are of the same size and very small compared to finite-difference grids used in conventional simulators. In practice, levels of mesh refinement in the continuum region of the simulator might be required, where the pore-scale region is coupled to the fine-scale continuum region, which is then coupled to a coarser continuum region. The simulator would also become very computationally difficult if extended to 3D. Computational efficiency was not a focus of this study, but future work will investigate and optimize 3D, time-dependent simulations in a parallel computing environment. Upscaling (a priori) to a conventional finite-difference grid is ideal and circumvents the need for hybrid simulation. However, we showed that significant inaccuracies can result because surrounding flow behavior, manifested in the form of boundary conditions, affects flow in the pore-scale models. One alternative is to use the a posteriori upscaled result. Although the multiscale simulation must be performed, it needs to be performed only once, and a permeability field is extracted for future simulations. Another option might involve solving each pore-scale model by coupling to a few surrounding subdomains (to provide boundary conditions), instead of the entire domain at once. This would be analogous to a local−global upscaling28−30 approach used in reservoir simulation. Suitable boundary conditions on a subset encompassing the subdomain might be estimated on the basis of the geographic location of the models with respect to the entire domain and its surrounding neighbors.
ACKNOWLEDGMENTS
This material is based upon work supported as part of the Center for Frontiers of Subsurface Energy Security, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0001114.
■
■
NOMENCLATURE K = permeability (L2) p = pressure (F/L2) q = flow rate (L3/T) r = radius (L) N = number of pores in the pore-scale model Ndof = number of degrees of freedom on mortar space α = coefficients of mortar basis functions ϕ = basis functions used in mortar space μ = viscosity (F L−2 T−1) REFERENCES
(1) Bryant, S.; Pallatt, N. Predicting formation factor and resistivity index in simple sandstones. J. Pet. Sci. Eng. 1996, 15 (2−4), 169−179. (2) Bakke, S.; Øren, P. E. 3-D pore-scale modelling of sandstones and flow simulations in the pore networks. SPE J. 1997, 2 (2), 136−149. (3) Valvatne, P. H.; Blunt, M. J. Water Resour. Res. 2004, 40, No. W07406. (4) Balhoff, M. T.; Wheeler, M. F. A predictive pore scale model of non-Darcy flow in porous media. SPE J. 2009, 14 (4), 579−587. (5) Joekar-Niasar, V.; Prodanović, M.; Wildenschild, D.; Hassanizadeh, S. M. Network model investigation of interfacial area, capillary pressure and saturation relationships in granular porous media. Water Resour. Res. 2010, 46, No. W06526. (6) Chen, S.; Doolen, G. Lattice Boltzman method for fluid flows. Annu. Rev. Fluid Mech. 1998, 30, 329−364. (7) Prodanovic, M.; Bryant, S. A level set method for determining critical curvatures for drainage and imbibitions. J. Colloid Interface Sci. 2006, 304, 442−458. (8) Fatt, I. The network model of porous media I. Capillary characteristics. Trans. Am. Inst. Min., Metall. Pet. Eng. 1956, 207, 144− 159. (9) Bryant, S. L.; Mellor, D. W.; Cade, C. A. Physically representative network models of transport in porous media. AIChE J. 1953, 39, 387−396. (10) Balhoff, M. T.; Thompson, K. E.; Hjortsø, M. Coupling pore networks to continuum models. Comput. Geosci. 2007, 33 (3), 393− 410. (11) Balhoff, M. T.; Thomas, S. G.; Wheeler, M. F. Mortar coupling and upscaling of pore scale models. Comput. Geosci. 2008, 12 (1), 15− 27. (12) Peterson, R. T.; Balhoff, M. T.; Bryant, S. Pore-scale modeling of the impact of surrounding flow behavior on multiphase flow properties. J. Multiscale Model. 2012, 3 (3), 109−131. (13) Sun, T.; Mehmani, Y.; Bhagmane, J.; Balhoff, M. T. Pore to continuum upscaling of permeability in heterogeneous porous media using mortars. Int. J. Oil, Gas Coal Technol. 2012, 5 (2/3), 249−266. (14) Scheibe, T. D.; Tartakovsky, A. M.; Tartakovsky, D. M.; Redden, G. D.; Meakin, P. Hybrid numerical methods for multiscale simulations of subsurface biogeochemical processes. J. Phys.: Conf. Ser. 2007, 78 (1), 1−5. (15) Tartakovsky, A. M.; Tartakovsky, D. M.; Scheibe, T. D.; Meakin, P. Hybrid simulations of reaction−diffusion systems in porous media. SIAM J. Sci. Comput. 2008, 30 (6), 2799−2816. (16) Koumoutsakos, P. Multiscale flow simulations using particles. Annu. Rev. Fluid Mech. 2005, 37, 457−487. (17) Chu, J.; Engquist, B.; Prodanovic, M.; Tsai, R. A multiscale method coupling network and continuum models in porous media I Single phase flow. Multiscale Model. Simul. 2012, 10, 515−549.
4. CONCLUSION A hybrid, multiscale sub-surface simulator has been developed that allows for direct substitution of pore-scale models for fine, continuum-scale grids. Here, a pore-scale region of over 7500 pore-scale network models was inserted around a well and coupled to an outer continuum scale region. Continuity of pressure and flux between subdomains was ensured using finiteelement mortars. Pore-level heterogeneity was captured, which would not be possible using a conventional, continuum simulator. Although network models were used in this work, other pore-scale approaches (e.g., Lattice Boltzmann) could be employed. A priori, direct upscaling can lead to inaccuracies for heterogeneous media in pressure and flux fields because the boundary conditions are arbitrary and do not reflect flow in surrounding media. However, an a posteriori, global approach showed that upscaled permeabilities were dependent upon the boundary conditions and resulted in a good approximation to the multiscale solution, both qualitatively and quantitatively. The simulator could be extended to model more nonlinear, dynamic flows that require pore-level modeling. Examples include the injection of carbon dioxide in sequestration strategies, non-Darcy flow, or matrix stimulation with acids. Future work should focus on these time-dependent applications and optimize computational efficiency in a parallel computing environment.
■
Article
AUTHOR INFORMATION
Corresponding Author
*Telephone: 512-471-3246. Fax: 512-471-4605. E-mail: balhoff@mail.utexas.edu. Notes
The authors declare no competing financial interest. 5835
dx.doi.org/10.1021/ef301003b | Energy Fuels 2012, 26, 5828−5836
Energy & Fuels
Article
(18) Arbogast, T.; Cowsar, L. C.; Wheeler, M. F.; Yotov, I. Mixed finite element methods on nonmatching multiblock grids. SIAM J. Numer. Anal. 2002, 37 (4), 1295−1315. (19) Wheeler, M. F.; Arbogast, T.; Bryant, S.; Eaton, J.; Lu, Q.; Peszynska, M.; Yotov, I. A parallel multiblock/multidomain approach for reservoir simulation. Proceedings of the Society of Petroleum Engineers (SPE) Reservoir Simulation Symposium; Houston, TX, Feb 14−17, 1999. (20) Peszynska, M.; Wheeler, M.; Yotov, I. Mortar upscaling for multiphase flow in porous media. Comput. Geosci. 2002, 6 (1), 73− 100. (21) Arbogast, T. Numerical Subgrid Upscaling of Two-Phase Flow in Porous Media; Texas Institute for Computational and Applied Mathematics (TICAM), University of Texas at Austin: Austin, TX, 1999; Technical Report 99-30. (22) Thomas, S. G.; Wheeler, M. F. Multiblock methods for coupled flow-transport and compositional flow through porous media Applications to simulation of transport of reactive species and carbon sequestration. Proceedings of the Society of Petroleum Engineers (SPE) Reservoir Simulation Symposium; The Woodlands, TX, Feb 21−23, 2011; SPE Paper 141824. (23) Al-Raoush, R.; Thompson, K. E.; Willson, C. S. Comparison of network generation techniques for unconsolidated porous media. Soil Sci. Soc. Am. J. 2003, 67 (6), 1687−1700. (24) Thompson, K. E.; Willson, C. S.; White, C. D.; Nyman, S.; Bhattacharya, J.; Reed, A. H. Application of a new grain-based reconstruction algorithm to microtomography images for quantitative characterization and flow modeling. Proceedings of the Society of Petroleum Engineers (SPE) Annual Technical Conference and Exhibition; Dallas, TX, Oct 9−12, 2005. (25) Jackson, M. D.; Valvatne, P. H.; Blunt, M. J. Prediction of wettability variation and its impact on flow using pore- to reservoirscale simulations. J. Pet. Sci. Eng. 2003, 39, 231−246. (26) Holden, L.; Nielsen, B. F. Global upscaling of permeability in heterogeneous reservoirs: The output least squares (OLS) method. Trans. Porous Media 2000, 40, 115−143. (27) Fredd, C. N.; Fogler, H. S. Influence of transport and reaction on wormhole formation in porous media. AIChE J. 1998, 44 (9), 1933−1949. (28) Chen, Y.; Durlofsky, L. J.; Gerritsen, M.; Wen, X. H. A coupled local-global upscaling approach for simulating flow in highly heterogeneous formations. Adv. Water Resour. 2003, 26 (10), 1041− 1060. (29) Chen, Y.; Durlofsky, L. J. Adaptive local-global upscaling for general flow scenarios in heterogeneous formations. Trans. Porous Media 2005, 62, 157−185. (30) Wen, X. H.; Durlofsky, L. J.; Chen, Y. Efficient threedimensional implementation of local-global upscaling for reservoir simulation. Proceedings of the Society of Petroleum Engineers (SPE) Reservoir Simulation Symposium; Houston, TX, Jan 31−Feb 2, 2005; SPE Paper 92965.
■
NOTE ADDED AFTER ASAP PUBLICATION This paper published August 13, 2012 without the Acknowledgments section. The correct version published August 20, 2012.
5836
dx.doi.org/10.1021/ef301003b | Energy Fuels 2012, 26, 5828−5836
