Modeling Three-Dimensional Anisotropic Structures of Reservoir Lithofacies Using Two-Dimensional Digital Outcrops

: Reservoir heterogeneity is a key geological problem that restricts oil and gas exploration and development of clastic rocks from the early to late stages. Existing reservoir heterogeneity modeling methods such as multiple-point geostatistics (MPS) can accurately model the two-dimensional anisotropic structures of reservoir lithofacies. However, three-dimensional training images are required to construct three-dimensional reservoir lithofacies anisotropic structures models, and the method to use reservoir heterogeneity model of fewer-dimensional to obtain a three-dimensional model has become a much-focused research topic. In this study, the outcrops of the second member of Qingshuihe Formation (K 1 q 2 ) in the northwestern margin of the Junggar Basin, which are lower cretaceous rocks, were the research target. The three-dimensional reservoir heterogeneity model of the K1q2 outcrop was established based on the unmanned aerial vehicle (UAV) digital outcrops model and MPS techniques, and the “sequential two-dimensional conditioning data” (s2Dcd) method was modiﬁed based on a sensitivity parameter analysis. Results of the parametric sensitivity analysis revealed that the isotropic multigrid simulations demonstrate poor performance because of the lack of three-dimensional training images, conditioning data that are horizontally discrete and vertically continuous, and distribution of lithofacies that are characterized by large horizontal continuities and small thicknesses. The reservoir lithofacies anisotropic structure reconstructions performed well with anisotropic multigrids. The simulation sequence of two-dimensional surfaces for generating the three-dimensional anisotropic structure of reservoir lithofacies models should be reasonably planned according to the actual geological data and limited hard data. In additional to this, the conditional probability density function of each two-dimensional training image should be fully utilized. The simulation results using only one two-dimensional section will have several types of noises, which is not consistent with the actual geological background. The anisotropic multigrid simulations and two-dimensional training image simulation sequence, proposed in this paper as “cross mesh, reﬁnement step by step”, e ﬀ ectively reduced the noise generated, made full use of the information from the two-dimensional training image, and reconstructed the three-dimensional reservoir lithofacies anisotropic structures models, thus conforming to the actual geological conditions. visualization, Y.Y.; supervision, Y.Y.; project administration, L.Z. and


Introduction
The spatial heterogeneity of reservoir parameters (porosity, permeability, and lithofacies) is called reservoir heterogeneity [1]. This heterogeneity affects the flow and accumulation of ground fluids (groundwater, oil, and gas) [2]. Reservoir heterogeneity modeling plays a crucial role in the exploration and development of underground oil and gas resources [3]. Traditional geological modeling methods can be divided into two categories: pixel-based, two-point geostatistics (such as truncated Gaussian simulations [4][5][6] and sequential indicator simulations [7][8][9]; and object-based methods [10][11][12][13], which encompass marked-point processes [14] and Boolean methods [15]). The former mainly uses the variogram of two points as a tool to reflect reservoir heterogeneity, which maybe faithful to the hard data, but cannot characterize the complex geometry (such as the shape of river channel) and spatial shape accurately. This method is more suitable for the marine strata [16][17][18][19][20][21][22][23]. The second method generates the spatial distribution of geometry shape, size, and direction of the simulated objected around the central point, according to the spatial law of the simulated objected [24]. Although in line with the experience of geologists, this method is not wholly faithful to hard data. The multipoint geostatistics method was proposed by Guardiano and Srivastava (1993) [25]. In this method, the probabilities are obtained directly from the training images, which does not have to agree with the hard data. The training image is a prior depiction geological model expressed the anisotropic structures of reservoir lithofacies and distribution pattern [3,26,27]. Based on the multipoint geostatistics of Guardiano and Srivastava (1993) [25], Strebelle (2000) [28] proposed the single normal equation simulation (SNESIM) algorithm, which uses a search tree to save the multipoint probability of scanned training images, and, based on this algorithm, applied the multipoint geostatistics effectively [3].
Many studies have researched the algorithms of multipoint geostatistics, such as SNESIM algorithm [28], SIMPAT algorithm [29], FILTERSIM algorithm [30], DISPAT algorithm [31]. Irrespective of the algorithm used, the simulated results are affected by the training image, especially the three-dimensional training image [32]. Incidentally, it is challenging to obtain a proper three-dimensional training image. The methods used to generate three-dimensional training images can be divided into three types: (1) Objected-based algorithm: directly used to build three-dimensional training images, the method is a numerical description of anisotropic structures of reservoir lithofacies [13,26,[33][34][35][36]. For example, fluvsim [13] and ellipsim [7] were applied by separate research groups to generate training images of the river channel and ellipsoid geometry, respectively. However, with this method, it is difficult to reconstruct the anisotropic structures of reservoir lithofacies and distribution patterns observed in the outcrops. (2) Process-based algorithm: used to construct three-dimensional training images [37], through analysis of background and actual geological data of the research area. The three-dimensional training image is constructed by the forward modeling techniques. For example, Zhang et al. (2017) used the Dionisos sedimentary forward modeling technique to establish a model of carbonate sedimentary rock [24]. Adopting this model as a training image and combining it with hard data and the SNESIM algorithm, a realistic simulation model was obtained. However, this method is rarely used in constructing a clastic reservoir lithofacies anisotropic structures model, and the process of multi-technology fusion is complex and difficult. (3) Three-dimensional training images based on two-dimensional training images and its statistical information [32,38,39]. Okabe and Blunt (2004) used this workflow to reconstruct the three-dimensional submicron-scale pore structure of carbonate, based on two-dimensional rock flakes, and improved the estimation of carbonate porosity [38]. Caers (2006) reconstructed the geological structures at different scales by using the statistical information of variogram and two-dimensional training images [39]. The "sequential two-dimensional conditioning data" (s2Dcd) workflow proposed by Comunian and Straubhaar (2012) [32] can be used to construct three-dimensional training images from two-dimensional training images, and it has been applied and promoted by several researchers [23,40]. The main idea behind this method is to use two-dimensional training images, which can represent the distribution considered by the researchers in the study area, to realize sequential simulation of two-dimensional slices. This sequence is defined by the user, and the conditioning data used are computed in the previous simulation steps and fed two-dimensional training images. These steps are repeated until the whole space is simulated [32]. Traditional geological methods [41] and the unmanned aerial vehicle (UAV) technology developed in recent decades [23,[42][43][44] can be used to accurately obtain two-dimensional digital outcrops and two-dimensional reservoir lithofacies anisotropic structures models. Three-dimensional training images are usually required to establish three-dimensional reservoir lithofacies anisotropic structures models of the outcrops through multipoint geostatistics, and the s2Dcd method solves this problem of using two-dimensional training images to construct three-dimensional reservoir lithofacies anisotropic structures models [32].
However, the application of the s2Dcd method still has some problems. For example, less noisy results were obtained by Pickel (2015) [23] using a single training image in the x-z path and transposing it along the y-z direction. Cannan (2016) [40] reduce noise of the results using the same method. Isotropic multigrid was used by Pickel (2015) [23], Cannan et al. (2016) [40] and Chen et al. (2018) [45] when utilizing or verifying this method. However, they found that, when using training images only in one direction and in isotropic multigrid, the multiple noise points and poor lithofacies continuity was caused in the simulation results. Moreover, the sensitivity of parameters used in the s2Dcd methods has not been explored much.
Therefore, this paper aims to measure through a case study the influence of sensitivity parameters and sequence for the sequential simulation of two-dimensional slices on the simulation results, and discuss the problems in the application of the s2Dcd method. This paper is organized as follows. First, Section 2 provides the geological background of the study area and basic knowledge on digital image acquisition. Section 3 consists of four parts: data collection and processing of UAV digital outcrops; the selection of multipoint geostatistics algorithm; architectural element analysis and classification; training image generation and multipoint geostatistics realization. Section 4 consists of multi-grid calculations and the sequence for the sequential simulation of two-dimensional surfaces, with the results interspersed throughout. Section 5 provides suggestions for future work. The method proposed by this paper will be very important to accurately predict the underground reservoir and estimate the reserves.

Geological Setting
The study area is located in the northwest border of the Junggar basin, the north of Wuerhe, Karamay city, Xinjiang provinces, China. The Qingshuihe Formations (lower Cretaceous), Carboniferous, and Quaternary were exposed in the study area (Figure 1a-c). The formation of the Junggar basin occurred during the Hercynian orogeny, its development being especially active in the late epoch of the Permian period (an epoch in the Hercynian orogeny). The Indosinian movement (257-205 Ma) in the late Triassic period formed the current tectonic pattern of the region. The Junggar basin in the Cretaceous was a gentle monoclinal structure inclined to the southeast and overlaid the Carboniferous strata to the north [46][47][48]. The Cretaceous strata, which are widely exposed in the study areas, can be divided four periods concerning sedimentary deposition: bottom conglomerate sedimentary period, river sedimentary period, shallow lake sedimentary period, and river and delta sedimentary period [46][47][48]. The object of this study is the second member of the lower Cretaceous Qingshuihe Formation (K1q2) (Figure 2), which is mainly composed of sandstone and conglomerate, interpreted to represent the deposits of a braided fluvial to deltaic system ( Figure 1d).

Digital Outcrop Model
In the past decade, the increasingly maturing UAV modeling technology has been widely used in scientific research [42,43,49]. In geological research, UAVs have several advantages in that they can be used to aerially study outcrops (such as a steep cliff or ravine) that cannot be studied manually due to the physical terrain conditions [42], measure and interpret the characteristics of outcrops of small-to medium-size sedimentary rocks [43], provide original data (digital outcrops) for the establishment of training images and two-dimensional reservoir lithofacies anisotropic structures models [23,44,49].
The workflow of generating a digital outcrop using a UAV can be broken down into the following steps: (1) obtain digital photographs with an overlap rate of over 60% of the outcrops, and record the location, elevation, and deflection angle in these photographs; (2) identify key points in each of the photographs; (3) match the key points in the photographs; (4) automatic aerial triangulation (AAT) and estimated camera pose based on beam-block adjustment (BBA); (5) process the directional photographs and obtain the cloud; and (6) generate the mesh and digital outcrops [43]. Several researchers have combined UAV modeling technology with multipoint statistical simulation method [20,23], and have proved that the three-dimensional surface model formed by UAV modeling can be used as the basic data for three-dimensional training image generation.

Digital Imagery Data Acquisition and Processing
The UAV, Phantom 4 Pro, developed by DJI, was used. The digital camera mounted on the cradle head of the UAV with a stable system is a one-inch CMOS camera with 20 million effective pixels. The camera scope rotates from −90 • to 30 • , with an angle jitter of ±0.02 • . The digital photographs of all Cretaceous outcrops in the study area were captured by the Phantom 4 Pro. In order to clarify the problems in the generation of three-dimensional models from two-dimensional training images, the outcrops marked with the red pentagram in Figure 1c were selected to provide the basic data for establishing of three-dimensional models in a later stage, and subsequent studies on the other outcrops will be published in the later. The first group digital photographs were captured by the UAV from a distance of approximately 2 m from the outcrop surface, to ensure that the small-scale sedimentary structure of the outcrop was clear. The second group digital photographs were captured from a distance of approximately 8 m from the outcrop surface, to ensure that the large-scale sedimentary structure of the outcrop and the connections between the three sections were clear. Multiple three-dimensional renderings of the outcrop were generated using the digital photographs on the Pix4D Mapper to create dense three-dimensional point clouds with overlying three-dimensional meshes. The digital photographs overlapped at the rate of 75%, the camera chooses the perspective of level (90 • ) to ensure the quality of the dense three-dimensional point clouds and three-dimensional meshes. Throughout the study area, the digital images were captured at an angle parallel to the outcrop surface to ensure that the true geometry of the outcrop was recorded. During the modeling, the uncertainty ellipses used to represent the image calibration effect in the Pix4D Mapper and the size of the uncertainty ellipses were same indicating that the digital images were calibrated properly and sufficient matches were found (more details can be found in these papers, e.g., Nieminski [42,43,49]). The number of point clouds generated on Pix4DMapper was 1,022,136. The accuracy of the model revealed by the three-dimensional rendering is 0.064 × 0.042 m. The boundary, proportion, length, and width of each lithofacies were tracked and measured by three-dimensional renderings in the reference coordinate system of Pix4DMapper. The UAV models are shown in Figure 3a.

Selecting Algorithms of Multipoint Geostatistics
Several multipoint geostatistics can be used, such as SNESIM algorithm [28], SIMPAT algorithm [29], FILTERSIM algorithm [30], DISPAT algorithm [31], and so on. However, these algorithms have their own advantages and disadvantages. SNESIM, for example, is the most typical of these algorithms, and is designed for modeling categories, e.g., facies distribution, and it is limited by the number of categories, which cannot exceed five. Moreover, it cannot be used for modeling a continuous property. SIMPAT [29] is a sequential simulation multipoint geostatistical algorithm with patterns. In SIMPAT, the idea of image reconstruction is applied in the field of geological modeling of reservoirs. For the estimation of unknown points, instead of the probability method, the similarity calculation method is used to compare the data modes in the training image and the data event at the point to be estimated. The problem affecting the spread of the SIMPAT algorithm is the calculated similarity between data events and all database styles. To solve these problems, Zhang et al. (2006) [30] proposed the FILTERSIM algorithms, which includes Filtersim_cate and Filtersim_cont. Filtersim_cate can be used for modeling categories, in which the number of the categories is more than five. Filtersim_cont can be used for continuous variable simulations, e.g., porosity and permeability of the reservoir. This type of algorithm requires linear filters, for different facies, and coding which is a challenging task. DISPAT algorithms [31] use multidimensional scaling analysis and k-means clustering to process the data style. The computational efficiency of DISPAT is much higher than SIMPAT and FILTERSIM; however, the multidimensional scaling analysis and k-means clustering in three-dimensional require large amounts of computational memory, which limits the practical application of DISPAT. In the application of the s2Dcd method [32], a multipoint geostatistics work platform for realizing SNESIM algorithms was required. The open source software, Stanford geostatistical modeling software (SGeMS) was used for this study. SGeMS can be downloaded from the website (http://sgems. sourceforge.net/). The classical algorithm (SNESIM) is used to establish the three-dimensional model or three-dimensional training image. There are two types of multiple grids simulation: isotropic and anisotropic in the SGeMS. In the previous studies, most researchers used 3th-or 4th-grade isotropic multiple grids [23,40,45]. In this paper, the isotropic and anisotropic multiple grids were all analyzed.

Rock Facies Analysis and Classification
The outcrops in Figure 3a were systematically sampled and cut into thin sections. The locations of these samples are shown in Figure 3b,d,e. The bottom of the outcrop is medium to coarse sandstone, which gradually changes to mudstone upward. The bottom of the outcrop is medium-scale trough cross-bedding, which gradually changes to high-angle cross-bedding upward and to low-angle cross-bedding. The top of the outcrops is mudstone with parallel bedding. There are five facies in the outcrops, the characteristics of the five facies are as follows (Figures 3 and 4): (1) Facies F1 is composed of coarse sandstone facies with high-angle cross-bedding and incised angle >15 • , interpreted to represent braided channel or braided bar of braided river. The microscopic characteristics of F1 are shown in W1 of Figure 4; the porosity of facies F1 is high (>15%), with rich bitumen. Samples W3, W4, W6, W8, W13, W16, W21, W23, W28, W30, and W34 belong to this category. (2) Facies F2 is composed of medium to fine sandstone facies with low-angle cross-bedding and an incised angle of <15 • . F2 is the margin of the braided river channel, and overlaps with the facies F1. The microscopic characteristics of F2 are shown in W2 of Figure 4 with rich plastic particles (>15%). The porosity of facies F2 is low (<6%). Samples W5, W15, W20, W20, W25, W31, and W32 belong to this category. (3) Facies F3 is composed of medium to fine sandstone facies with low-angle cross-bedding and macroscopic characteristics similar to facies F2. The microscopic characteristics of F3 are shown in W13 of Figure 4 with rich gypsum cement (>15%). The porosity of facies F3 is low (<3%). Samples W14, W19, and W27 belong to this category. (4) Facies F4 is composed of medium to fine sandstone facies with small-scale through cross-bedding.
The thickness of this type of facies is relatively thin, which may be due to the late deposition of braided river. The microscopic characteristics of F4 are shown in W22 of Figure 4 with rich calcium cement (>15%). The porosity of facies F4 is low (<3%). Samples W9, W10, W11, W18, W24, W35, and W37 belong to this category. (5) Facies F5 is composed of medium thick mudstone facies and may have interchannel mud deposits.
The microscopic characteristics this sample are shown in W7 of Figure 4, without bituminous soil. The porosity of facies F5 is low (<1%). Samples W29, W36, and W38 belong to this category.
The outcrop (Figure 3a) was studied according to the field methodology proposed by Miall (1988) [41], and this information was used as the basis for classification. Figure 3a shows the digital outcrops and tracking interface of facies, Figure 3b-e is the section interpretation results corresponding to A-A', B-B', and C-C' in Figure 3a. The lengths of the three corresponding sections are 837, 603, and 603 cm, respectively. The height is 312 cm. Facies F2 and F1 include the 4th order bounding surfaces corresponding to the macroform growth increment. Therefore, lithofacies of F1, F2, F3, F4, and F5 are the separated 5th order bounding surfaces corresponding to the channel. According to the results of reservoir architectural analysis, distribution of the five lithofacies in the three sections is described. In order to model the anisotropic structures of reservoir lithofacies through multipoint geostatistics, the surfaces which bounding the five lithofacies in the UAV digital model were vertically projected onto the corresponding plane. Sections A-A' and C-C' are defined as sections in the x direction, and section B-B' is defined as a section in the z direction. The surfaces of UAV digital model of x and z direction were constructed, and these three surfaces were the basic data for multipoint geostatistics simulation (Figure 5a).

Parameters and MPS Realization
Based on the number of lithofacies, the number of categories was defined as five, where each category corresponds to one lithofacies. The target proportion of the five lithofacies in the simulated section was defined the same as the five lithofacies in the corresponding parallel two-dimensional training image. For example, the target proportion of the five lithofacies in the simulated section, which is the parallel z section, was defined as: F1:F2:F3:F4:F5 = 0.35:0.15:0.21:0.07:0.22. The servo system factor was set to 0.5. To improve the simulation quality, the nodes in the search template was set to 115. The conditional data are the data points on the intersection lines between the known sections or the pre-simulated sections and the simulated sections. The search ellipsoid should be designed prior to running the simulation of multipoint geostatistics on SGeMS. Six parameters are included while designing the search ellipsoid: three search radii, r max , r med , and r min and three rotation angles, α, β, and γ. Before rotation, r max is the search radius in the y-axis, r med is the search radius in the x-axis, and r min is the search radius in the z-axis. The three rotation angles, α, β, and γ in the search ellipsoid are rotated about z', x', and y', which are the axes of the ellipsoid.
The lithofacies distribution in the outcrops shows that the lithofacies extension direction is close to x-axis or z-axis, and the angle is approximately 1 • . Therefore, if the simulated section is parallel to the x-y section, then α = 89 • , and β = γ = 0. If the simulated section is parallel to the z-y section, then β = 89 • , α = γ = 0. In addition, when the previous sections were simulated, the conditioning data were horizontally discrete and vertically continuous. The search ellipsoid should include all data on the sides, so the r max was designed to be 603, while the r med was 200, the r min was 190.
The UAV-aided models upon vertical projection and digitization can provide three lateral sections as the conditioning data and two-dimensional training image (Figure 5a). Three training images were constructed with 1 × 1 cm grids. There are 837 grids in the x-axis, 603 grids in the z-axis, and 312 grids in the y-axis. Before generation of three-dimensional training image, a two-dimensional simulation of length 837 on section C-C' was simulated, with section A-A' as the training image, and section C-C' as the conditioning data, since section C-C' and section A-A' have different lengths. The result is shown in Figure 5b. At this point, the basic data for generating the three-dimensional training images or reservoir lithofacies anisotropic structures models were fully prepared. The five lithofacies varied greatly along the y-axis. For example, the top and bottom of the model were composed of mudstones, while the middle of the model was composed of five facies. Moreover, all the two-dimensional training images were in z and x directions. Therefore, the y-coordinate was used as an auxiliary variable. In this way, a three-dimensional training image or three-dimensional reservoir lithofacies anisotropic structures models was generated, with a series of two-dimensional simulation in the x and z directions.

Multiple Grids Simulation and Its Influence on Simulation Results
In order to accurately capture the large-scale lithofacies anisotropic structures, the concept of multiple grids simulation was proposed by Tran (1994) [50]. Section 2 from the model simulated without multiple grids (Figure 6a) shows that the lithofacies anisotropic structures in the x-coordinate or z-coordinate could not be captured accurately, while those in the y-coordinate could be captured only fairly (Figure 6b). Therefore, multiple grids were incorporated to improve the quality of the generation three-dimensional models. With the isotropic multiple grids, the simulation results remained unchanged when the grade was less than 2 (g < 2). The results with g = 2, 3, 4, and 5 are shown in Figure 6c. When g was greater than 2, the larger-scale lithofacies anisotropic structures in the y-coordinate were captured much better. However, the anisotropic structures in the z-coordinate with the same isotropic multiple grids for the same values of g were captured poorly. The actual data in the training image (section z) reveals that the length of the z-coordinate is approximately 1.9 times that of the y-coordinate. The conditioning data was on the side of section (2), and the connectivity of lithofacies in the z-coordinate was stronger than that in y-coordinate. The isotropic multiple grids cannot capture larger-scale anisotropic structures in the z-coordinate and y-coordinate simultaneously. Therefore, because of the anisotropic connectivity of the lithofacies and distribution of the conditioning data, the isotropic type is not feasible for the analysis. Moreover, the anisotropic multiple grids are constructed to improve the quality of simulation results. For section (2), the number of anisotropic multiple grids in the x-coordinate by design is 1, because the number of grids in the x-coordinate is 1. Based on the results of previous simulation, it can be seen that the larger-scale lithofacies anisotropic structures in the y-coordinate were captured well when the coarse grid was greater than 2. The grade of anisotropic multiple grids in the y-coordinate by design was 2, that in the coarse grid was 2, and that in the fine grid was 1. After a sensitivity analysis of the anisotropic parameters, the grade of the anisotropic multiple grids in the z-coordinate were designed as 2, the coarse as 21, and the fine grid as 1. The results using anisotropic multiple grids are shown in Figure 6d-f. When the anisotropic multiple grids were used, larger-scale structures were captured well in the z-coordinate and y-coordinate, and the connectivity of lithofacies in the z-coordinate and y-coordinate was strong. However, the results using anisotropic multiple grids are still not realistic.
Firstly, there is still a lot of noise in the simulation results, especially in the middle position of Figure 6d.
Secondly, anisotropic multiple grids make low utilization of the conditioning data. The conditioning data of section (2) is on two intersecting lines (where sections A-A' and section C-C' intersect with section (2)). The conditioning data on the intersection between section C-C' and section (2) were well-utilized in all anisotropic multiple grids. Thus, it can be concluded that section (2) reproduced the lithofacies continuity and larger-scale lithofacies anisotropic structures better on the right sides (z = 0) (Figure 6f). However, the conditioning data on the intersection between section A-A' and section (2) were not well-utilized in any of the anisotropic multiple grids, which means that the continuity of lithofacies could not be reproduced on the left side (z = 602) (Figure 6e). The reason was that the expansion factor in the z-coordinate was 21 (note that "21" is the expansion factor of the multiple grids, not the 21st multiple grid), and the number of grids in the z-coordinate was 603, the first number of z-coordinate was z = 0, and the label of the left-most coarse grid is z = 588. However, the conditioning data on the intersection between section A-A' and section (2) was on the label of z = 602. This set of conditioning data appeared to have been discarded during the simulation, so the continuity of lithofacies could not be reproduced (Figure 6e), and it does not play any substantial role in the reconstruction of lithofacies anisotropic structures. For section (2), there were only two sets of conditioning data, which is not acceptable. To solve this problem, it was necessary to gain a better understanding of the transferring and discarding of conditioning data in the multiple grid simulation, which will be described in detail in the next section.

Data Relocation in Multiple Grids
The analysis in the previous sections has proved that larger-scale lithofacies anisotropic structures can be reconstructed using 2nd grade anisotropic multiple grids. However, the utilization of hard data is still insufficient. Figure 7a,c shows the utilization of conditioning data when using 2nd anisotropic multiple grids. After the z-coordinate is roughened with the expansion factor of 21, the coarse grid is closest to the conditioning data on the line where z = 602 is z = 588. The distance between the conditioning data and the coarse grid is 15, and the data events with conditioning data on the sides of the section have a low repetition rate. Therefore, the conditioning data on the intersection between section A-A' and section (2) was discarded in the coarse grid simulation. This set of conditioning data is used only in the fine grid simulation. The conditioning data on the intersection between section A-A' and section (2) did not play a controlling role in the reconstruction of the larger-scale lithofacies anisotropic structures. Moreover, this resulted in the lithofacies in Figure 6e being discontinuous. In order to solve this discontinuity in the simulation process, the conditioning data on the z = 602 was translated to z = 588 (Figure 7b), to improve the control effect of this set of conditioning data on the simulation results. The utilization of conditioning data after relocation are shown in Figure 7b,d. The conditioning data were fully used, and the noise in the simulation results was suppressed efficiently (Figure 8a,b). The results shown in Figure 8a,b reveal that the larger-scale lithofacies anisotropic structures on the z-coordinate and y-coordinate, and the connectivity of facies F1, which is the oil-bearing layer at the bottom of the model were well reconstructed. The lithofacies continuity at the connection between sections A-A' and C-C' was also well reconstructed. However, it was difficult to understand whether the conditioning data translated were consistent with the actual geological conditions. The authors reviewed Tran's (1994) [50] definition of multiple grids. When the grid is coarsened, the conditioning data is often hard to locate. Thus, the conditioning data on the coarse grid should be moved to the nearest coarse grid. This method can effectively improve the simulation results when the conditioning data are relatively discrete and evenly distributed. The conditioning data in section (2) was not discrete or evenly distributed; thus, manual relocation of the conditioning data was required.
However, this method still has some disadvantages. When the conditioning data are used to translate 15 grids, the real distance of the translated grids is only 15 cm, as the size of each grid is 1 × 1 cm. Hence, the method is fairly reliable for this grid size. However, when the size of each grid is 1 × 1 km, the real distance of the translated grids is 15 km, which is often not allowed for the actual geological conditions.

Sequence of Simulation Surfaces and Its Influence on Simulation Results
When Comunian and Straubhaar (2012) [32] proposed the s2dcd method, the constraint on the simulation sequence of two-dimensional training image was defined as follows: "If the dimensions of the simulation domain (hereinafter is the length of the x and z coordinates) are comparable, then a uniform cross grid sequence of simulation can be used; If x y or z y, the sequence of simulation can be in a single direction, while the training image in the other direction is mainly used as the conditioning data." The results achieved by Comunian and Straubhaar (2012) [32] indicate that the diversity of sequences of simulation surfaces has a great impact on the simulation results. In order to reduce the noise, Pickel et al. (2015) [23] incorporated a single training image and transposed it along both two-dimensional paths. Although Pickel et al. (2015) [23] using the method were able to reduce the noise in the three-dimensional models, the simulation results show that their method cannot reflect the heterogeneity in different directions due to river migration. Due to these reasons, the training image can only be used when it is parallel to the simulation surface. The actual data reveal that the length of the x-coordinate is approximately 1.4 times that of z-coordinate. To study the influences of simulation surface sequence on the results, this section compared the results of using a unidirectional training image and bidirectional training image with different simulation surface sequences. This paper designs four simulation surfaces sequences, which is shown in Figure 9. Figure 9a illustrates the two-dimensional simulations were 835 times in the z direction with the z directional training images only. The simulated results are shown in Figure 10c,d. Figure 9b shows two-dimensional simulations were executed 601 times in the x-y direction with the x-y directional training images only. These results are shown in Figure 10a,b. Figure 9c shows that, firstly, two-dimensional simulations were executed 39 times on the coarse grids in the z direction with the z directional training images; secondly, two-dimensional simulations were executed 28 times on the coarse grids in the x-y direction with the x directional training images; and thirdly, two-dimensional simulations were executed 796 times on the fine grids in the z direction with the z directional training images. The simulated results are shown in Figure 10e,f. The sequence of simulation surfaces as shown in Figure 9d was proposed in this study. Regardless of the relationship between the number of grids in the x-coordinate and z-coordinate, the simulated three-dimension space can be divided into several secondary spaces by several surfaces in the z direction and x direction. Then, the secondary spaces can be divided into four equal thirds by several surfaces in the z direction and x direction, and the thirds can be divided into four equal fourth spaces by several surfaces in the z direction and x direction, and so on. The sequence of simulation surface is also the sequence of construction the secondary space. To generate the conditional data for the next simulation, which are located in the coarse grids, the label of the surface which was used to construct the secondary space showed a multiple of the maximum expansion factor. The advantage of this method is that the conditioning data in the next simulation will not be offset, and the phenomenon of lithofacies discontinuity shown in Figure 6a can be avoided. For example, the number of grids in the x-coordinate was 837, and the number of grids in the z-coordinate was 603, and the maximum expansion factor in the z-coordinate or x-coordinate was 21; thus the three-dimensional space can be divided in to nine secondary spaces by five surfaces (the five surfaces are as follows: surface on which x = 210; surface on which x = 420; surface on which x = 630; surface on which z = 210; surface on which z = 420). Moreover, every secondary space can be divided into four equal third spaces by two surfaces. Consequently, the smallest space was 1 grid (x-coordinate) × 1 grid (y-coordinate). The sequence was named, "cross mesh, refinement step by step." The simulated results with this sequence are shown in Figure 10e,f. The following observations can be made from the simulation results:

1.
When using a single directional training image (x-y direction or z-y direction), the output results showed that the lithofacies continuity and lithofacies anisotropic structure were well-reconstructed in the training image direction, but they were poorly reconstructed in the vertical direction of training image. For example, when the image in the z-y direction was used as the training image, and the training image in the x-y direction was used as the conditioning data, the lithofacies continuity and anisotropic structure in the z-y direction were well-reconstructed, but poorly in the two vertical directions. These results are shown in Figure 10c,d. Similarly, this was observed when the training image in the x-y direction was used as the training image and the training image in the z-y direction was used as the conditioning data. The results are not in accordance with the actual situation (Figure 10a,b). These two sets of simulations confirmed the direction in which the training image was used, and that the pattern of lithofacies in the training image direction was well-reconstructed, while the reconstruction in the other two directions was poor. Therefore, when the three-dimensional reservoir lithofacies anisotropic structures models were well-constructed with a unidirectional training image, the patterns in the other two directions could be barely reproduced.

2.
When using the bidirectional sections (x and z directions) as the training image, the output results exhibited a significant improvement in the reconstruction of lithofacies heterogeneous structures. Figure 10e,f shows the simulation results using a bidirectional training image, and its simulated sequence has been described above. This simulation sequence significantly improves the heterogeneous structure of the lithofacies in a direction parallel to the z direction when compared to those shown in Figure 10c,d, using a unidirectional training image. However, the phenomenon of lithofacies discontinuity still appears in the middle of the coarse grid in the z direction. In addition, lithofacies continuity and lithofacies anisotropic structure in the z direction were reconstructed poorly (Figure 10e,f). Noise in the simulation results persisted after the latter were improved, especially in the z direction. Although this simulation sequence uses bidirectional training images, the pattern of lithofacies in the x direction was just reconstructed in the section of the coarse grids. The pattern of lithofacies of fine grids between coarse grids in the z direction could not be reconstructed.
In view of this drawback, the sequence of simulation surface shown in Figure 9d was proposed to improve the lithofacies continuity in all directions in the simulation results. The results shown in Figure 10g,h, using the "cross mesh, refinement step by step" show well-reconstructed lithofacies continuity and lithofacies anisotropic structure in three directions. The results also show that the conditional probability density function and conditioning data in the three training images were fully utilized. The reconstruction of lithofacies distribution patterns in the x-z direction is better than the other three sequences. Moreover, the spatial distribution of lithofacies F1 reveals that the braided river flows mainly in the north-south direction, and that oil and gas mainly occur in the main braided channel or the braided bar of the braided river channel.
Therefore, the workflow of the s2Dcd method in the application process was modified as follows: 1.
Based on the parameter sensitivity analysis of the two-dimensional simulation, the general parameters and expansion factors of anisotropic multiple grids in all directions were determined.

2.
Based on the expansion factor, the distribution of conditioning data, and the lengths of the three axes, the sequence of simulation surface was determined. 3.
The determined sequence of simulation surface was executed. 4.
The simulation surface was selected.

5.
Conditioning data in the simulating surface was chosen. 6.
Steps c-f were executed in a loop. 8.
The method was terminated, and the output simulation results were obtained.
However, there are still some problems in establishing the three-dimensional reservoir lithofacies anisotropic structures models using two-dimensional outcrops or training images although the simulated sequence proposed in this paper reduces the noise in the simulation results, and makes the simulation results close to practical geological scenarios, the noise cannot be completely eliminated in the simulation results. In addition, the verification of three-dimensional models of the anisotropic structures could not be performed. Most research works have attempted modify the algorithms and simulation methods for generating three-dimensional models from two-dimensional training images. The methods are validated by comparing the nuclear magnetic resonance (NMR) results with the simulation results in the porous media. The spatial distribution of porous media has low categories (rock skeleton and pore space) and is relatively uniform in three directions (x-y direction, z-y direction, and x-z direction). These are different from the structures in the outcrops, which are usually different in the three directions (e.g., along the direction of flow of the river, vertical to the direction of flow of the river, and parallel to the direction of flow of the river), and there are more than two types of lithofacies.
In addition, the influence of three-dimensional models of anisotropic structures of reservoir lithofacies on the fluid flow and accumulation has not been discussed yet. Although this working fluid flow can be used to establish a three-dimensional reservoir lithofacies anisotropic structures model, the multipoint algorithm SNESIM used in this paper cannot simulate the continuous phase. The simulation of fluid flow and accumulation in this model could be a challenge to overcome in the future. These problems should be solved in subsequent research work.