## Abstract

Small-scale (mm to m) sedimentary structures (e.g. ripple lamination, cross-bedding) have received a great deal of attention in sedimentary geology. The influence of depositional heterogeneity on subsurface fluid flow is now widely recognized, but incorporating these features in physically-rational bedform models at various scales remains problematic. The current investigation expands the capability of an existing set of open-source codes, allowing generation of high-resolution 3D bedform architecture models. The implemented modifications enable the generation of 3D digital models consisting of laminae and matrix (binary field) with characteristic depositional architecture. The binary model is then populated with petrophysical properties using a textural approach for additional analysis such as statistical characterization, property upscaling, and single and multiphase fluid flow simulation. One example binary model with corresponding threshold capillary pressure field and the scripts used to generate them are provided, but the approach can be used to generate dozens of previously documented common facies models and a variety of property assignments. An application using the example model is presented simulating buoyant fluid (CO_{2}) migration and resulting saturation distribution.

## Introduction

Two longstanding challenges in flow and transport problems have been the accurate representation of geologic heterogeneity in numerical models and quantifying the influence of that heterogeneity on fluid flow^{1,2}. Here heterogeneity refers to mm- to m-scale variability in grain sizes and architecture (spatial organization) associated with clastic depositional processes as manifested in descriptive structures. There is a long history of study of depositional sedimentary structures at these scales and their influence on fluid flow has been known for decades^{3,4,5,6,7}. Critical aspects of geological heterogeneity related to small-scale depositional features have been documented for a variety of considerations^{8,9,10,11,12,13}. These types of studies highlight the importance of understanding the influence of depositional heterogeneity at relatively small scales and the challenges involved in adequately representing such heterogeneity for various applications.

The digital models presented here are a variant of traditional geostatistical representations of depositional facies^{1,14}, which many geoscientists agree were a vast improvement over homogeneous/isotropic representations and quite advantageous, but often fall short of effectively representing the small scale details in their most recognizable and descriptive forms. The method for generating bedform architecture models (BAM) presented here is most closely associated with geometric approaches^{15,16} in that the models are meant to be geologically plausible, depositionally realistic, and visually satisfying from a geologic perspective. The focus is on relating the depositional processes and resulting characteristic descriptive architecture to the distribution of material properties (i.e. grain size, permeability, porosity, threshold capillary pressure). Geostatistical techniques are often used for populating large-scale models for applications such as reservoir flow simulation^{17} and hydrogeological applications^{18}. While geostatistical techniques continue to evolve in complexity^{19}, their ability to generate depositionally realistic small-scale representations of standard geologic features remains a challenge. The approach here differs from the sedimentary process forward model approach that requires complex computational fluid dynamics, including turbulence^{20}. Other attempts at representing sedimentary architecture at small scale include utilization of natural specimens or outcrops in two and three dimensions^{21,22,23,24}.

The approach presented here extends a well-documented and widely utilized set of codes originally presented by Rubin^{25} and later published by Rubin and Carter^{26} as a set of MATLAB routines. David Rubin pioneered the digital representation of clastic depositional features, linking our understanding of sediment transport processes, bedform morphology, and descriptive sedimentology. The ability to generate bedform images and visualize bedform deposition and migration represented a major advancement in sedimentary geology and has been utilized by sedimentologists and stratigraphers to understand sedimentary processes and to interpret structures observed in the field and subsurface cores^{27}. The original Rubin^{25} computer program was based on a geometric model, using sine curves to generate surfaces that mimic well-documented bedforms. The numerical modification of up to three sine curves simulated the migration and erosion of bedforms, defined by 75 user-definable geometric parameters^{26}. While this approach does not mimic hydrodynamic processes related to deposition, it does produce readily recognizable features (i.e. similar to outcrop and modern environments) that can be related to depositional processes. The original output consists of three-dimensional block diagram images representing vertical and horizontal sections of the bedforms (Fig. 1) along with polar plots of bedform migration directions. Rubin and Carter^{26} subsequently revised and rewrote the code in MATLAB and published scripts that, in addition to the block diagrams, generate outputs with 3D bedform topography as well as animation sequences that illustrate bedform deposition over time. As illustrated in Table 1, bedforms can be grouped based on descriptive aspects such as planform shape (2D, 3D), behavior through time (variable, invariable) and crest orientation (transverse, oblique, longitudinal). Even though 3D illustrations (digital images) could still be generated with the updated codes, they do not produce 3D geocellular volumes. This presented an opportunity to modify the animation sequence protocol to produce 3D digital models which can then be further manipulated for various computational applications. This provides yet another breakthrough related to Rubin’s formative work, by allowing broad application of realistic 3D sedimentologic geocellular models tied to well-documented geologic understanding for exploring a wide range of traditional research topics in sedimentary geology and fluid flow simulation.

## Results and Discussion

One advantage of characterizing small scale heterogeneity through the set of codes derived from Rubin’s work is that it encourages users to understand how input parameters affect the final model structure. This presents an opportunity to more closely tie descriptive understanding of clastic sediments and their depositional environments with physical properties and fluid flow performance.

Geostatistical techniques are unlikely to be replaced by the methods presented here, but rather they can be further integrated into the models presented. For example, the binary model output representing the depositional laminae can be used as a training image or conditioning fabric for other stochastic approaches, such as multiple point statistics (MPS)^{28}. Such enhancements are likely to produce even more realistic models than those generated here, potentially including normal or inverse grading of bed sediment.

We have not yet explored a variety of model modifications that might be pursued. For example, the digital model populated with petrophysical parameters (e.g. permeability) could be smoothed to provide a more continuous property distribution (as opposed to random sampling from a PDF). This would provide more locally correlated values that are probably more representative of the continuous and smooth changes in properties that are observed in natural specimens^{23}. In addition, the incorporation of stochastic noise for laminae position assignment could be employed to modify the strictly periodic representation of current models that result from using sine and cosine functions for generating bedform architectures. Such modifications are likely to present even more geologically-sound representations of depositional facies.

Modification of the Rubin and Carter^{26} codes allows for the generation of realistic digital models of 3D depositional bedform architecture. These models combine a deterministic bedform architecture component mimicking realistic crossbedding geometries with stochastic variability of petrophysical properties. Another advantage of this approach is that it allows consideration of domain sizes larger than the core plugs typically used for laboratory flow experiments, where small sizes may not fully capture depositional architecture. An example of one binary model (BAM #36) with corresponding *P*_{
th
} (threshold capillary pressure) field is provided with the script used to generate it, but the approach is applicable for dozens of previously documented BAMs (and input parameter files) in Rubin^{25}. One example application illustrates the use of three different BAMs representing ripple-laminated fluvial sediment for simulating buoyant CO2 migration through *P*_{
th
} field with increasing textural contrast between matrix and laminae, resulting in volumetric fractions of trapped CO_{2} ranging between 2.7% and 69.5% of the model domain.

### Example model

The scripts (m-files) provided enable the user to generate binary digital geocellular models of any of the 79 figures from Rubin and Carter^{26} and to populate each model with any of the facies from Beard and Weyl^{29} (Fig. 2). An example binary model (provided as a MATLAB workspace variable) with a corresponding threshold capillary pressure field that can be generated using the scripts provided is available in the supplementary material (*vert_xsxn_all.mat*). The BAM presented is # 36 of Rubin and Carter^{26}, and the block diagram image that is part of the output is shown in Fig. 1. With assigned cell dimensions of 2 × 2 × 2 mm^{3}, the model is 90 × 60 × 35 cm^{3}. The model matrix is populated with moderately sorted upper coarse sand (MUCSa) and the laminae with moderately sorted upper fine sand (MUFSa), and threshold capillary pressure values are assigned as described above.

### Practical implication

To provide an example of how the datasets generated can be used, an application is presented illustrating the use of the example model for simulating buoyant fluid migration of supercritical CO_{2} dominated by capillary forces. The movement of carbon dioxide in the subsurface is an important consideration for projects that inject CO_{2} into subsurface reservoirs^{31,32}. Far from injection wells, and after injection ceases, buoyancy and capillary forces are known to control flow behaviour, while small scale heterogeneities strongly influence saturation distributions^{33}. As a matter of fact, small-scale depositional features have been shown to have significant influence on CO_{2} migration both numerically^{5,34,35,36} and experimentally^{37,38}. Due to the capillary-dominated flow assumptions, the numerical simulations can be carried out in a matter of seconds using the invasion percolation algorithm^{39,40}. We use Permedia software to model buoyant migration of CO_{2} from a planar source at the bottom of the model until the CO_{2} plume spans the entire model (i.e., at percolation). The goal of these simulations is to understand the effect of *P*_{
th
} contrast between matrix and laminae as well as the geometry of the crossbedding on the volumetric fraction of trapped CO_{2} plume of each BAM. Figure 3 shows simulation results for three different textural contrasts applied to three different BAMs (#5, #36, and #63). Results from 200 equally-probable realizations of the *P*_{
th
} field (using the same architecture, but different realizations of property assignment) indicate that volumetric fractions of trapped CO_{2} plume in these models is expected to range between 2.7% and 69.5% of the model domain. This high variability provides a realistic expectation for CO_{2} saturations in this type of clastic material for the subsurface conditions considered. This topic is more fully and rigorously explored with more bedform architecture models and facies assignments in Trevisan, *et al*.^{41}.

## Methods

### Description of the modified code

The starting point for code development is the MATLAB code by Rubin and Carter^{26}. This set of scripts can be used to create images of depositional features with known morphology and interpreted genesis, and can be tied to specific geologic examples that are well documented in the literature^{16,25,27}. In the modified scripts provided here, a 3D geocellular model is intuitively built from adjacent 2D vertical sections parallel to a vertical face represented in the block diagrams produced. This is needed because the original script does not explicitly have as part of the output a workspace variable in MATLAB with values representing 3D assignment of matrix and laminae cells.

The primary modification of the original *NewDunes.m* function is a new outer loop that allows a 3D dataset to be generated from adjacent 2D sections (essentially serial slicing vertical planes) of the standard 3D block diagram output. For advanced users who want to understand the procedure which allows a 3D model output, the *for* loop beginning at line 32 of the modified code (*Model_Generator.m*) incrementally changes the input variable PHASEF, which is described by Rubin and Carter^{26} as the bedform phase (in degrees). Changing PHASEF essentially specifies the horizontal position of the bedform within the fixed reference frame of the block diagram edges. By systematically changing this variable and retaining the vertical face of the block diagram output for each loop, serial sections of the model can then be concatenated to build a 3D digital model.

Beyond the addition of simple looping, the modified code takes advantage of a built-in MATLAB function that is perhaps non-intuitive for 3D geologic model building. The modified script uses the ‘vec2mtx’ (line 106 of *Model_Generator.m*) function, which is used in mapping applications and is available through MATLAB’s Mapping Toolbox. More specifically, the ‘vec2mtx’ function converts X-Y vectors (polygons) to a specified regular data grid of binary values (zero or one). Practically, it allows for a continuous representation of depositional laminae in a gridded format, similar to pixelated (raster) representation of a map of political boundaries. This function is applied to each vertical section of the model during the looping of PHASEF described above. This mapping function is needed because the workspace variable representing the laminae on the block diagram faces is initially stored as a three column array representing coordinates of the laminae locations, and is not a regularized grid. For additional information about ‘vec2mtx’, the reader is referred to the MATLAB documentation and numerous online examples.

For users familiar with the routines by Rubin and Carter^{26}, the script published here contains comments on each line that represents a modification of the original script or new entry. The new script will produce an image of the model, but all scripting related to visual animation of bedforms in the original script has been eliminated for clarity and compactness. Since the original script contains standard comment notation, new and/or modified lines are identified in the new scripts by the use of four comment characters (%%%%). The comments provide a description of the purpose of the modification or new entry made in each line. As such, the described modifications are fairly minor.

Model domain size (the number of cells in the orthogonal horizontal directions) can be modified by changing the ‘grid density’ variable (line 28 of *Model_Generator.m*) and by defining the area of interest (lines 95, 103 and 105 of *Model_Generator.m*). Higher values of grid density represent the same architecture at increasing resolution, and model sizes of millions of cells are not uncommon for low values of grid density (e.g. 5). The size of the output 3D models is typically assigned considering the wavelengths of the input variables, as defined in the variable input file (e.g. SPCNGF; line 3 of any *Fig##.m* parameter file from Rubin and Carter^{26}). This input parameter is defined as the bedform wavelength in the first set, and similar input variables correspond to the second or third bedform set (SPCNGS, SPCNGT, where the final letter corresponds to the set number). Maximum model domain size capability for the modified code was not fully explored, or a variety of other aspects of model generation, as our ultimate goal was standard facies model utilization for flow simulation. The model domain sizes created by the default variables in the figure parameter files of Rubin and Carter^{26} were sufficient for our intended purposes. Other potential (unexplored) modifications are included in the Discussion section.

The output 3D model in the modified scripts does not have prescribed cell physical dimensions. Model dimensions are user-specified after generating the model, and cell dimensions need to be defined by the user in subsequent model applications depending on the knowledge of sedimentary structures, which can be informed by Rubin^{25}, examples in literature^{42}, and personal experience. Obviously there is subjectivity in cell dimension assignment, but some assignments would produce models with features that are not readily found in natural depositional settings at those scales. That is, sedimentary structures formed by known depositional processes have typical physical dimensions, which should be reflected by reasonable assignment of cell size (and resultant model domain dimensions), based on the user’s knowledge and experience. In the model application presented here, cubic sub-volumes of approximately one million cells with cell dimensions of 0.202 × 0.202 × 0.202 m^{3} are extracted from the 23,625,000-cell model produced by *Model_Generator.m*. The procedure to generate these representative elementary volumes (REV) is explained in appendix A of Trevisan *et al*.^{41}.

The output of *Model_Generator.m* includes an image of the model domain in the vertical Y-Z plane (Fig. 4A) including the defined area of interest for the model, a fence diagram of the model showing three orthogonal planes (Fig. 4B), as well as a MATLAB workspace variable (“vert_xsxn_all”) with lamina cells represented by ones and matrix cells represented by zeros.

Subsequently, the entirely new script *Value_Assign.m* needs to be run to assign petrophysical properties such as grain sizes and threshold capillary pressure (*P*_{
th
}) to the individual model cells (e.g. MATLAB workspace variable named “vert_xsxn_all_trunc_Pth”; Fig. 4C). The assignment of any non-binary physical properties is user defined and can be extended beyond the methods presented here by editing the code, but in the method described here is fundamentally related to specific geologic facies as described below. That method can be adopted, or another preferred method could be employed by the user.

### Property assignment

In addition to generating a synthetic 3D model on a Cartesian grid, the set of codes populates the binary model with petrophysical properties such as grain size, permeability, or threshold capillary pressure. In order to do this, matrix and lamina cells of the binary models are assigned with textural classes characterized by a range of grain sizes and grain sorting.

Beard and Weyl^{29} classified 48 clastic facies of artificially mixed unconsolidated sands into eight sand grain-size subclasses and six sorting groups, whereas Fig. 2 also includes a ninth column with silt. Following previous work of Meckel, *et al*.^{43} we assign the matrix and lamina cells of the output 3D model with the lithofacies shown in Fig. 2. With the purpose of representing realistic depositional environments, grain sizes assigned to matrix facies are always coarser than grain sizes assigned to the laminae. Assignment of identical lithofacies to both laminae and matrix will create a rather homogeneous model, whereas varying amounts of textural contrast between matrix and laminae will regulate the baffle effect exerted by the laminae. For the example application presented, the matrix is considered to be made up of moderately sorted upper coarse sand (MUCSa), while the laminae are composed of, from low to high contrast, moderately sorted upper medium sand (MUMSa), upper fine sand (MUFSa), and upper very fine sand (MUVFSa). The assignment of facies to matrix and laminae is user-defined by screen prompt during execution of the script *Value_Assign.m*. The options for these assignments range from median grain sizes of upper coarse sand to upper coarse silt, and from extremely well sorted to very poorly sorted.

After matrix and lamina facies have been assigned to the binary model, grain size values for each cell are sampled from probability density functions (PDF) corresponding to each lithofacies. The sampling from the probability distribution of grain sizes and corresponding petrophysical properties occurs in the *Value_Assign.m* script via built-in random number generator and is a necessary step to account for the high uncertainty of natural geological heterogeneity. The mean and standard deviation of each PDF are related to the median grain diameter and sorting category, respectively. Each of the facies types shown in Fig. 2 has lognormally distributed grain diameters^{44} and is described by the median diameter (*d*_{50}) and the Trask sorting coefficient (*S*_{
o
}). The Trask sorting coefficient is defined as the square root of the ratio of the 25^{th} to the 75^{th} percentile of the cumulative density function of the lognormal distribution:

where *d* is grain diameter.

Since the grain sizes are often log normally distributed^{44}, the logarithm of the grain size (*D*) will be normally distributed and is given as:

Then, substituting 2 into 1, the Trask sorting coefficient can be re-written as:

For a normal distribution, the quantile function is defined as:

where *µ* is the mean, *σ* is the standard deviation, *erf* is the error function, and *p* is the cumulative probability.

The mean of the normal distribution is given by the logarithm of the median of the lognormal distribution:

Substituting 4 and 5 into 3, we solve for the standard deviation *σ* of the normal distribution. Then, using the calculated mean and standard deviation, we calculate the inverse CDF for cumulative probabilities from 0 to 0.999, which gives the values of the normal distribution. Since these are logarithmic values, the antilog gives back the lognormally distributed grain diameters. These values are then substituted in equation 6 to obtain lognormally distributed *P*_{
th
} values. The lognormal distribution of grain diameters is thus reverse constructed only by using the median and standard deviation values of each of the facies reported in the original paper by Beard and Weyl^{29}.

However, if used without any geological knowledge, the results can assign the model cells unnaturally (and irrationally, when considering associated depositional flows) large or small grain size values (outliers in the long tails of lognormal distributions) with subsequent impact on other calculated petrophysical properties (e.g. permeability). So careful scrutiny of model properties is necessary while assigning values from the lognormal distribution.

Beard and Weyl^{29} prepared facies samples by combining sand grains of different sizes such that they were lognormally distributed. But they report that the sand grains used ranged from 0.4 mm to 9 mm in diameter. The same size limits are thus here enforced while populating the models with the codes presented. This is implemented in the program by truncating the distribution between those two grain diameters whenever necessary. The truncation comes into effect only when the distribution is very wide (poorly and very poorly sorted). The code has a provision to change this as per user’s requirements, if there is a preferred method for assigning such properties.

Any grain size distribution obtained in the previous section can be used to calculate equivalent petrophysical parameters like porosity, permeability and capillary entry pressures. For the specific application presented below (i.e. capillary-dominated fluid migration) we populate the model with *P*_{
th
} values.

The capillary entry pressure *P*_{
th
} (in kPa) is obtained from the grain size using equation 6 from Berg^{45}:

where 16.3 is a geometric constant, *γ* is interfacial tension in N/m, and *d*_{50} is median grain diameter, in millimeters. For the example simulation provided we consider 0.03 N/m for CO_{2}-H_{2}O system, representative for reservoirs at approximately 10 MPa and 35 °C^{46} (~1.5 km depth). These properties are relevant for geological carbon sequestration.

### Data availability statement

All the required MATLAB scripts and one example output dataset generated during the current study are included in this published article (and its Supplementary Information files).

## Change history

### 06 March 2018

A correction to this article has been published and is linked from the HTML and PDF versions of this paper. The error has been fixed in the paper.

## References

- 1.
Koltermann, C. E. & Gorelick, S. M. Heterogeneity in sedimentary deposits: A review of structure-imitating, process-imitating, and descriptive approaches.

*Water Resour Res***32**, 2617–2658, doi:10.1029/96wr00025 (1996). - 2.
Soltanian, M. R. & Ritzi, R. W. A new method for analysis of variance of the hydraulic and reactive attributes of aquifers as linked to hierarchical and multiscaled sedimentary architecture.

*Water Resour Res***50**, 9766–9776 (2014). - 3.
Cushman, J. H.

*Dynamics of fluids in hierarchical porous media*. (Academic Press, 1990). - 4.
Pryshlak, T. T., Sawyer, A. H., Stonedahl, S. H. & Soltanian, M. R. Multiscale hyporheic exchange through strongly heterogeneous sediments.

*Water Resour Res***51**, 9127–9140 (2015). - 5.
Gershenzon, N. I.

*et al*. Influence of small-scale fluvial architecture on CO_{2}trapping processes in deep brine reservoirs.*Water Resour Res***51**, 8240–8256, doi:10.1002/2015wr017638 (2015). - 6.
Soltanian, M. R., Ritzi, R. W., Huang, C. C. & Dai, Z. Relating reactive solute transport to hierarchical and multiscale sedimentary architecture in a Lagrangian‐based transport model: 1. Time‐dependent effective retardation factor.

*Water Resour Res***51**, 1586–1600 (2015). - 7.
Soltanian, M. R., Ritzi, R. W., Huang, C. C. & Dai, Z. Relating reactive solute transport to hierarchical and multiscale sedimentary architecture in a Lagrangian‐based transport model: 2. Particle displacement variance.

*Water Resour Res***51**, 1601–1618 (2015). - 8.
Dawe, R. A., Caruana, A. & Grattoni, C. A. Immiscible Displacement in Cross-Bedded Heterogeneous Porous Media.

*Transp. Porous Media***87**, 335–353, doi:10.1007/s11242-010-9687-4 (2011). - 9.
Pickup, G. E., Ringrose, P. S., Jensen, J. L. & Sorbie, K. S. Permeability tensors for sedimentary structures.

*Math Geol***26**, 227–250, doi:10.1007/bf02082765 (1994). - 10.
Weber, K. J. Influence of common sedimentary structures on fluid flow in reservoir models.

*J Petrol Technol***34**, 665–672 (1982). - 11.
Luo, X. R.

*et al*. Effects of carrier bed heterogeneity on hydrocarbon migration.*Mar Petrol Geol***68**, 120–131, doi:10.1016/j.marpetgeo.2015.08.015 (2015). - 12.
Kuo, C. W. & Benson, S. M. Numerical and analytical study of effects of small scale heterogeneity on CO

_{2}/brine multiphase flow system in horizontal corefloods.*Adv Water Resour***79**, 1–17, doi:10.1016/j.advwatres.2015.01.012 (2015). - 13.
Ringrose, P. S., Sorbie, K. S., Corbett, P. W. M. & Jensen, J. L. Immiscible flow behaviour in laminated and cross-bedded sandstones.

*J Petrol Sci Eng***9**, 103–124, doi:10.1016/0920-4105(93)90071-L (1993). - 14.
Haldorsen, H. H. & Damsleth, E. Stochastic Modeling.

*J Petrol Technol***42**, 404–412 (1990). doi:SPE-20321-PA. - 15.
Ramanathan, R.

*et al*. Simulating the heterogeneity in braided channel belt deposits: 1. A geometric‐based methodology and code.*Water Resour Res***46**(2010). - 16.
Scheibe, T. D. & Freyberg, D. L. Use of sedimentological information for geometric simulation of natural porous media structure.

*Water Resour Res***31**, 3259–3270 (1995). - 17.
Deutsch, C. V. & Tran, T. T. FLUVSIM: a program for object-based stochastic modeling of fluvial depositional systems.

*Comput Geosci-Uk***28**, 525–535, doi:10.1016/S0098-3004(01)00075-9 (2002). - 18.
Koch, J. & Nowak, W. Predicting DNAPL mass discharge and contaminated site longevity probabilities: Conceptual model and high‐resolution stochastic simulation.

*Water Resour Res***51**, 806–831 (2015). - 19.
Mustapha, H., Dimitrakopoulos, R. & Chatterjee, S. Geologic heterogeneity representation using high‐order spatial cumulants for subsurface flow and transport simulations.

*Water Resour Res***47**(2011). - 20.
Olsen, N. R. B. Three-dimensional CFD modeling of self-forming meandering channel.

*Journal of Hydraulic Engineering***129**, 366–372 (2003). - 21.
Basumallick, S. Size differentiation in a cross-stratified unit.

*Sedimentology***6**, 35–68 (1966). - 22.
Bayer, P., Comunian, A., Hoyng, D. & Mariethoz, G. High resolution multi-facies realizations of sedimentary reservoir and aquifer analogs.

*Sci Data***2**, 150033, doi:10.1038/sdata.2015.33 (2015). - 23.
Meckel, T. A. Digital Rendering of Sedimentary-Relief Peels: Implications for Clastic Facies Characterization and Fluid Flow.

*J Sediment Res***83**, 495–501, doi:10.2110/jsr.2013.43 (2013). - 24.
Comunian, A., Renard, P. & Straubhaar, J. 3D multiple-point statistics simulation using 2D training images.

*Comput Geosci-Uk***40**, 49–65, doi:10.1016/j.cageo.2011.07.009 (2012). - 25.
Rubin, D. M. In

*Society of Economic Paleontologists and Mineralolgists*Vol. 1 187 (Concepts in Sedimentology and Paleontology, 1987). - 26.
Rubin, D. M. & Carter, C. L. Bedforms and Cross-Bedding in Animation. (U.S. Geological Survey, 2005).

- 27.
Kocurek, G. Interpretation of Ancient Eolian Sand Dunes.

*Annual Review of Earth and Planetary Sciences***19**, 43–75, doi:10.1146/annurev.earth.19.1.43 (1991). - 28.
Strebelle, S. Conditional simulation of complex geological structures using multiple-point statistics.

*Math Geol***34**, 1–21, doi:10.1023/A:1014009426274 (2002). - 29.
Beard, D. C. & Weyl, P. K. Influence of Texture on Porosity and Permeability of Unconsolidated Sand.

*Am Assoc Petr Geol B***57**, 349–369 (1973). - 30.
Folk, R. L. & Ward, W. C. Brazos River Bar: a study in the significance of grain size parameters.

*Journal of sedimentary petrology***27**, 3–26 (1957). - 31.
Soltanian, M. R.

*et al*. Simulating the Cranfield geological carbon sequestration project with high-resolution static models and an accurate equation of state.*International Journal of Greenhouse Gas Control***54**, 282–296 (2016). - 32.
Soltanian, M. R., Amooie, M. A., Dai, Z., Cole, D. & Moortgat, J. Critical Dynamics of Gravito-Convective Mixing in Geological Carbon Sequestration.

*Sci Rep***6**, 35921, doi:10.1038/srep35921 (2016). - 33.
Cavanagh, A. J. & Haszeldine, R. S. The Sleipner storage site: Capillary flow modeling of a layered CO

_{2}plume requires fractured shale barriers within the Utsira Formation.*International Journal of Greenhouse Gas Control***21**, 101–112, doi:10.1016/j.ijggc.2013.11.017 (2014). - 34.
Green, C. P. & Ennis-King, J. Effect of Vertical Heterogeneity on Long-Term Migration of CO

_{2}in Saline Formations.*Transp. Porous Media***82**, 31–47, doi:10.1007/s11242-009-9498-7 (2010). - 35.
Li, B. X. & Benson, S. M. Influence of small-scale heterogeneity on upward CO

_{2}plume migration in storage aquifers.*Adv Water Resour***83**, 389–404, doi:10.1016/j.advwatres.2015.07.010 (2015). - 36.
Soltanian, M. R.

*et al*. Dissolution Trapping of Carbon Dioxide in Heterogeneous Aquifers.*Environ Sci Technol*(2017). - 37.
Krevor, S. C. M., Pini, R., Li, B. X. & Benson, S. M. Capillary heterogeneity trapping of CO

_{2}in a sandstone rock at reservoir conditions.*Geophys Res Lett***38**, doi:10.1029/2011gl048239 (2011). - 38.
Trevisan, L.

*et al*. Experimental analysis of spatial correlation effects on capillary trapping of supercritical CO_{2}at the intermediate laboratory scale in heterogeneous porous media.*Water Resour Res***51**, 8791–8805, doi:10.1002/2015wr017440 (2015). - 39.
Wilkinson, D. & Willemsen, J. F. Invasion Percolation - a New Form of Percolation Theory.

*J Phys a-Math Gen***16**, 3365–3376, doi:10.1088/0305-4470/16/14/028 (1983). - 40.
Carruthers, D. In

*Multidimensional Basin Modeling*Vol. 7 (ed S. Duppenbecker and R. Marzi) 21–37 (AAPG/Datapages Discovery Series, 2003). - 41.
Trevisan, L., Krishnamurthy, P. G. & Meckel, T. A. Impact of 3D capillary heterogeneity and bedform architecture at the sub-meter scale on CO

_{2}saturation for buoyant flow in clastic aquifers.*International Journal of Greenhouse Gas Control***56**, 237–249, doi:10.1016/j.ijggc.2016.12.001 (2017). - 42.
Boggs, S.

*Principles of sedimentology and stratigraphy*. 5th edn, (Pearson Prentice Hall, 2012). - 43.
Meckel, T. A., Bryant, S. L. & Ganesh, P. R. Characterization and prediction of CO

_{2}saturation resulting from modeling buoyant fluid migration in 2D heterogeneous geologic fabrics.*International Journal of Greenhouse Gas Control***34**, 85–96, doi:10.1016/j.ijggc.2014.12.010 (2015). - 44.
Folk, R. L. A Review of Grain-Size Parameters.

*Sedimentology***6**, 73–93, doi:10.1111/j.1365-3091.1966.tb01572.x (1966). - 45.
Berg, R. R. Capillary Pressures in Stratigraphic Traps.

*Aapg Bulletin-American Association of Petroleum Geologists***59**, 939–956 (1975). - 46.
Chiquet, P., Daridon, J. L., Broseta, D. & Thibeau, S. CO

_{2}/water interfacial tensions under pressure and temperature conditions of CO_{2}geological storage.*Energ Convers Manage***48**, 736–744, doi:10.1016/j.enconman.2006.09.011 (2007).

## Acknowledgements

This work could not have been undertaken without the prior seminal work of Rubin^{25} and Rubin and Carter^{26}. The work presented here was supported entirely by the Center for Frontiers of Subsurface Energy Security, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award # DE-SC0001114.

## Author information

### Affiliations

### Contributions

T.A.M. designed research and generated the modified MATLAB code for model generation and P.G.K. generated the MATLAB code for property assignment. All authors contributed to the research results, written main manuscript text, and figures. All authors reviewed the manuscript for content and accuracy.

### Corresponding author

## Ethics declarations

### Competing Interests

The authors declare that they have no competing interests.

## Additional information

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

A correction to this article is available online at https://doi.org/10.1038/s41598-018-21903-y.

## Electronic supplementary material

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Meckel, T.A., Trevisan, L. & Krishnamurthy, P.G. A method to generate small-scale, high-resolution sedimentary bedform architecture models representing realistic geologic facies.
*Sci Rep* **7, **9238 (2017). https://doi.org/10.1038/s41598-017-09065-9

Received:

Accepted:

Published:

## Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.