A New Modeling Method for the Sandbody Architecture of Braided River Reservoirs: A Case Study from HG Formation, N Gas Field

Reservoir architecture refers to the geometry, sizes, directions, and stacking patterns between the formation units and the interbeds of different levels of reservoirs [1, 2]. Based on the method and technique of Miall’s river facies architecture analysis, the fine characterization of the geometry size of different grade reservoir constituent units and their interrelationships within the pan-connected giant thick braided river complex sand body is an effective means to reveal the reservoir’s internal seepage barrier distribution patterns and improve the field development effect.

At present, scholars have achieved a lot of results and understanding in braided river reservoir sand types and characteristics [3–5], quantitative relationships of scale parameters [6–8], and sand stacking patterns [9, 10] through modern sedimentation, field outcrops, and flume experiments. Modern sedimentation and simulation experiments show that the top layer elevation of the sand body deposited in the channel bar within the braided river is higher than that of the braided river channel [11, 12]. Chen et al. [13] measured the size of several modern sandy braided river channel bars in Yukon River, Jamuna River, and Rakaia River. The length of the bars ranged from 300 to 1000 m, with an average of 800 m. The width of the bars ranged from 100 to 600 m, with an average of 280 m. Kelly [14] used 22 modern braided rivers (or flume experimental dates) and 34 ancient outcrop data to establish the equations between the width of the channel bar in sandy braided river (⁠$wb$⁠) and the depth of single channel (⁠$dm$⁠) and the length of single channel bar (⁠$lb$⁠) and channel bar width (⁠$wb$⁠): $wb=11.413dm1.4182$⁠; $lb=4.9517wb0.9676$⁠. With the development of braided river reservoir exploration [15–17], how to quantitatively characterize results of detailed braided river reservoir architecture analysis in the 3D geological models, so as to achieve the purpose of characterizing the internal inhomogeneity, has gradually become the focus and difficulty of braided river reservoir characterization. Braided river oil and gas reservoirs are widely developed in Cenozoic terrestrial basins. As the most important conformational unit of the sandy braided river reservoir, the 3D quantitative simulation of its interface and conformational unit, which is accurately characterized in the model, can improve the prediction accuracy of the remaining oil distribution and better guide the later adjustment and tapping of the field. Scholars at home and abroad have conducted a lot of exploration and research on this. Some practical and feasible technical methods have been developed [18–23] to promote the application of practical oilfield development.

Yao et al. [24], Yin et al. [25], and Bai [26] used object-based modeling method with progressive constraints layer by layer to establish the internal conformational model and parametric model of the braided river. Object-based modeling is able to define the geometry of each sandbody. However, there are still shortcomings in characterizing the architecture patterns and conditionality of the sandbodies in the braided river reservoir. Liu et al. [27] and Xu et al. [28] used the established training images as the basis to establish the braided river sedimentary microphase model by multipoint geostatistics. This method can characterize the interrelationship between braided channels and channel bars. However, its simulation results are discrete, so it cannot better portray the scale and morphology of the local plutons. Sun et al. [29] and Liu et al. [30] used sequential indicator simulation with combined human-machine postprocessing method for braided river 3D architecture simulation. The simulation process of this method is controlled by several parameters, such as the variance of the variance function and azimuthal angle. However, the geological body’s spatial geometry and distribution pattern are not well simulated. Therefore, this paper proposes a new modeling method for braided river reservoirs. A series of architecture pattern, morphological parameters, quantitative relationships, and well data obtained from a priori geological understanding is fully applied to quantitatively portray the geometry and spatial distribution of the channel bars. The comparison of modeling effects in the actual work area shows that the new modeling method simulates a 3D reservoir architecture model that is more consistent with the actual geological understanding and improves the accuracy of the model.

The conceptual model of the geologist’s architecture is the basis for defining the basic elements of the simulation. Take the braided river sandbodies architecture as an example. The braided river has a stable river channel. In the wide river valley, the sandbodies distribute continuously. The main sediments are channel-fill deposits and channel bar sandbodies. A braided river depositional architecture pattern with repeated superposition of several phases of channel bars and braided river channel in planar and vertical directions [26] (Figure 1). Field outcrop and sink experiments have shown [10, 31] that there are 6 major categories and 21 sand contact styles in the braided river reservoir (Figure 2). The lateral overlay reflects the lateral migration of the river channel. Single sand bars represent the accretion of primary source sand dams and the composite accretion of primary and secondary sources. The composite sand bar reflects the lateral migration and growth of the channel bar. Lateral splicing reflects the lateral splicing and intersection of primary and secondary source sandbodies in space. The different microfacies scales and vertical stacking characteristics within the braided river channel can be used as an important reference for modeling the subsurface braided river reservoir architecture. The establishment of a suitable a prior architecture model for the braided river sandbody is a prerequisite for the architecture modeling of the channel bar.

The basic process for the sandbody architecture of braided river reservoirs proposed in this paper is as follows: (1) first, without distinguishing the braided channel and the channel bar, the two are combined into sandstone, and a sand mudstone model is established; (2) starting from a known well, according to the thickness of the channel bar delineated on the well, combined with the thickness distribution of a single channel bar, and determine whether it is a multiphase channel bar. Then, expand outward according to the size and geometry of the channel bar. Completion of the simulation of the channel bar on the well; (3) determining whether the simulated channel bar reaches the set ratio. If not, the grid is randomly selected between wells. Based on the thickness of the vertical sandbody at that location, determine if a channel bar needs to be simulated. If so, simulate and update the channel bar percentage. The process is repeated until the given percentage of channel bar is reached.

During the simulation, only the center location, thickness, and aspect ratio parameters of each channel bar, and the spatial topological relationships between individual channel bars, between channel bar and the channel need to be adjusted. These adjustments are made within the constraints of the architecture patterns. According to the sedimentological principle, the lateral width, thickness, and length of the sedimentary sandbody are governed by the hydrodynamic conditions at the time of deposition. There are certain quantitative statistical relationships between the width and thickness and the length and width of the sandbody. The size of the sandbody can be determined based on the geometric parameters of the sandbody of the channel bar. The thickness of the channel bar is expressed as the vertical radius $c$ of the ellipsoid. The width-thickness ratio of the channel bar is expressed as the ratio of the short axis diameter 2a to c of the ellipsoidal plane. The length-to-width ratio of the channel bar is expressed as the ratio of the long and short diameters of the ellipsoid plane, 2b to 2a. These data can be obtained from actual work area information as well as from a priori knowledge. The parameters related to the simulation of this method are shown in Table 1.

The specific steps of the algorithm are as follows:

The adjustment principle of the channel bar of Step 7 is as follows.

When the channel bar does not satisfy the conditions given in Step 7, the first consideration is that the orientation of the channel bar does not match the local river orientation. Therefore, the channel bar is rotated by a certain angle to find the most appropriate direction of the channel bar. If the condition is still not satisfied after finding the direction, the coordinates of the center point of the channel bar are then shifted within a certain range based on the most appropriate direction. After a certain position adjustment, the channel bar will be more in line with the actual geological model.

The specific adjustment methods are divided into Option 1 and Option 2:

• (1)

  Option 1. Prefer strategy 1. Center point coordinates remain unchanged. Rotating the channel bar within the given range to find the angle that best matches the current river segment orientation.

  • (a)

    Dividing the rotation angle interval by interpolation method. Based on the initially set angle range of the channel bar, the rotation angle range is widened. The upper limit increases angle₁, and the lower limit decreases angle₂. Then, the new range is divided into $n$ angle intervals by setting the interpolation value angle₃

  • (b)

    Determining the rotation angle of the channel bar. Randomly select an angle interval. Then taking a random angle θ within this interval as the rotation angle of the channel bar. Deleting the angle interval that has already been selected

  • (c)

    Check the result. Return to Step 7 to check if this channel bar satisfies the given condition. If so, go to Step 8 or Step 10, otherwise return to (b). If none of the $n$ intervals satisfy the condition, then perform (d).

  • (d)

    Screening the angle that best matches the current channel segment. If none of the $n$ intervals meet the condition, the direction with the smallest proportion of nonsand body facies in the channel bar of the $n+1$ channel bar simulation results (the above $n$ angles +1 default source direction simulation result) is set as the direction of this channel bar. Option 2 is then executed

• (2)

  Option 2. Choosing Option 2 when Option 1 is not satisfied. Moving the center point coordinates in the given range.

  • (a)

    Determine the range of movement. Moving the coordinates of the center point of the channel bar in a circle with the center point p (⁠$x$⁠, $y$⁠, and $z1$⁠) on the well as the center point of the channel bar and a radius r₁

  • (b)

    Calculate the scalar distances from all grid points (Figure 5, p₁, p₂, and p_n) in the circle to the center point p of the bar at known well points. The distances are sorted from smallest to largest

  • (c)

    Determining the coordinates of the center point of the channel bar. The grid point with smaller scalar distance is preferred as the center point

  • (d)

    Checking the results. If the condition is satisfied, Step 8 or Step 10 is executed, otherwise return to (b). If the condition is not satisfied after traversing all the grids in the circle, then (e) is executed

  • (e)

    If the condition is not satisfied after traversing all the grids in the circle. Then, setting the width of the bar equal to min_ac$∗$c and the length equal to min_ba$∗$a. In the simulation results of all grids in the circle as the center point, the point with the smallest proportion of nonsand bodies in a single bar is used as the center point of the moved bar. Then, perform Step 8 or Step 10

The unconditional braided river sand petrographic model was established using an object-based method. At this point, no distinction is made between braided river flow channel and channel bar. Based on this model, a new modeling method for the sandbody architecture proposed in this paper is used to establish an architecture model containing the channel bar. The feasibility of the new modeling method is tested.

On the basis of braided river sandstone phase model, the grid coordinates (⁠$x$ and $y$⁠) are chosen randomly in the plane. The relationship between the thickness of the continuous sand body on the grid column determined by this grid and the thickness threshold of the channel bar is used to divide the number of superimpositions of the channel bar in the vertical direction. Figures 6 and 7 show the simulations of single-phase and multiphase superposition of the channel bar in the vertical direction. In this case, the single braided channel is 1000 m wide and 8 m thick.

Simulation of the single-phase channel bar in the vertical direction (Figure 6), the thickness of the continuous sand body at the simulated point in the vertical direction satisfies the condition of simulating only one channel bar. Different simulation realizations were obtained by varying different values of the thickness and width-to-thickness ratio parameters of the channel bars. For example, (1) if the thickness of the Heartland Dam is 7-8 m and the width-to-thickness ratio is 60 : 1 to 62 : 1, the simulated realization is obtained (Figure 6(a)). (2) If the thickness of the heart beach dam is 5 ~ 6 m and the width-to-thickness ratio is 50 : 1 ~ 60 : 1, a simulated realization is obtained (Figure 6(b)). (3) If the width-to-thickness ratio is further reduced, the channel width is sufficient to accommodate multiple channel bars, and the simulated realization is obtained (Figure 6(c)). It can be seen that if the thickness and width of the braided channel sand body are fixed. For larger values of the thickness parameter of the channel bars, only one period of channel bars is simulated in the vertical direction. The smaller parameter value of the width-to-thickness ratio, the greater the number of channel bars that can be accommodated in the lateral direction.

Figure 7 shows the superimposed simulation of a multiphase heartland dam. The superimposed migration of the multiperiod river channel results in a large thickness of continuous sand at the simulated points, which satisfies the conditions for simulating multiperiod channel bars in the vertical direction. The size and shape of the channel bars are controlled by setting the values of the following three main parameters. For example, multiple simulations based on the same braided channel sandstone phase model with the new algorithm yield multiple realizations. In this case, the main parameters of the channel bars are 3 ~ 8 m in thickness, 20 : 1 ~ 50 : 1 in width-to-thickness ratio, and 2 : 1 ~ 4 : 1 in length-to-width ratio. The braided channel thickness value is more than two times the lower limit of the thickness of the channel bars. Therefore, multi-phase river superposition can be simulated at suitable locations. When the thickness of the continuous sand body of the river at the simulated point is fixed, the thickness of the channel bar is the main parameter affecting the number of superimposed periods of the channel bar at that position. The smaller value, the more the number of periods of the channel bar in the vertical direction. Based on the J-directional slices achieved by the simulation in Figure 7, it can be seen that the channel bars simulated by this algorithm have the same pattern of stacking in the vertical direction as the braided river sand body configuration distribution.

The channel bars simulated based only on parameters such as thickness, center point location, width-to-thickness ratio, and aspect ratio of the channel bars may have the following deficiencies: (1) some or even most of the channel bars are located in mudstone; (2) there is a superposition relationship between the channel bars simulated by different wells or grid columns, and the independence of the channel bars cannot be maintained; (3) the direction of the channel bar does not match the direction of the current river segment, such as the channel bar in the circle of Figure 8(a). Therefore, this paper proposes 2 strategies to adjust the location of the channel bar to correct the above existing problems. First, set the parameter p1, and the smaller the p₁, the lower proportion of nonsandstone phases within a single channel bar (Figures 8(a)–8(c)). It is judged whether the proportion of nonsand body phases within the simulated channel bars has reached p₁. If not, rotate the channel bar by strategy one, or combine strategy one and strategy two to rotate and then translate the channel bar so that the simulated channel bar meets the given conditions. If not, rotate the channel bar by strategy one, or combine strategy one and strategy two to rotate and then translate the channel bar so that the simulated channel bar meets the given conditions.

Rotating channel bar. It is not necessary to change the quantitative parameters of the geometry of the channel bar and the coordinates of the center point. It is only necessary to rotate the channel bar by a certain angle and find the direction (angle value) of the channel bar that best fits the current river section. The angle interval of the channel bar is expanded by setting the lower limit decrease value angle₁ and the upper limit increase value angle₂ of the angle interval that can be rotated. That can satisfy the change of the angle of the channel bar in the local river section. The smaller interpolation angle₃ within the angle rotation interval, the more angle rotation intervals are available for random selection of the channel bar. Figure 8 shows an example where the angle interval of the channel bars is set to 20° and 45°, and the local change of channel bar’s direction due to the river bend may be beyond 20° and 45°. By setting different angle₁ and angle₂, the angle interval of the channel bars are expanded to (20°- angle₁ and 45° + angle₂), and the channel bars are rotated within this range. As can be seen in Figure 8, angle₃ remains unchanged, and the orientation of the channel bars matches more closely with the direction of the river as the range of angles over which the channel bars can be rotated is expanded. Meanwhile, this strategy can reduce the proportion of nonsandstone phase in the channel bars relatively quickly for channels with large local curvature (Figure 9).

Translating channel bar. When strategy one does not yet satisfy the given conditions, just move the channel bar in a circle of radius r₁ on the basis of the most suitable direction determined by strategy one. During the movement, the well is always located inside the channel bar. For example, on both sides of the river near the bank, the thickness of the sand body decreases toward the bank, and the proportion of nonsand body phase within the simulated channel bar is large. The requirement cannot be achieved by rotation alone, and at this point, the channel bar is shifted in the direction of the river’s centerline.

From the above, it can be seen that parameters such as thickness, width-thickness ratio, and aspect ratio of the channel bar can determine its size. And parameters such as p₁, r₁, and angular rotation range determine the spatial location of the channel bar. By adjusting these parameters, the distribution pattern of the internal architecture of the braided river reservoir sand body in the longitudinal direction can be better characterized, and the quantitative portrayal of the geometry and scale of the channel bar can be realized.

The West Lake depression is located in the central part of the eastern depression of the East China Sea shelf basin and is the largest Cenozoic oil and gas bearing depression in the East China Sea region [33]. It is adjacent to the Diaoyu Island uplift to the east, the Sea Reef uplift to the west, and the Fujiang Depression and the Diaobei Depression to the north and south, respectively. It is generally northeast trending. It is divided into five tectonic units from west to east, including the western slope zone, the western subconcave, the central inversion tectonic zone, the eastern subconcave, and the eastern fault zone (Figure 10) [34, 35].

The N gas field is located in the central inversion zone of the Xihu Depression [36]. According to the well data [37], the HG Formation (Figure 10) in the study area develops a giant thick sand body (the main producing layer). The target formation is located in the central H₃b section, which is subdivided into 5 sublayers. In this paper, the H₃b₁ sublayer is selected as an example for testing and comparing the new method. Zhu et al. have combined a large number of cores, logs, and seismic data with a comprehensive analysis of the regional depositional background (paleodepth, paleotopography, and source conditions), sediment flow mechanism, and sand body spreading characteristics. It is concluded that H₃b₁ mainly develops braided river delta deposits. The flow line of the submarine diversion channel is mainly in the northeast-southwest direction. It is supplied by dual sources in the northeast and east directions, of which the northeast direction is the main source direction [38]. The wells are all located near the channel. The occurrence of multiphase divergent channel superposition leads to a continuous spreading of giant thick sand bodies along the submerged divergent channels [39, 40]. The sandy content is comparatively channel bars that are more developed, followed by muddy deposits and river stagnation deposits. Different sedimentary bodies cut, spliced, and staggered increase the spatial anisotropy and nonhomogeneity. The main sedimentary microphase channels and channel bars are frequently migrated and superimposed, forming a series of composite channels and composite channel bar. This increases the spatial complexity and nonhomogeneity of the reservoir. Channel bars are the dominant phase zones developed in the sweet spot, and it is crucial to accurately characterize the channel bars in the lithofacies model.

The thickness of the heart beach dam in this area ranges from 0.6 m to 11 m, and the width ranges from 70 m to 400 m. The width-to-thickness ratio of the central beach mainly ranges from 25 to 62. The channel bars aspect ratio is concentrated between 2 : 1 and 4 : 1, with an average value of about 3 : 1.

In order to finely characterize the development of the channel bars, the average step length of the grid is set to $50m×50m$ in the plane, and the step length of the grid is set to 0.6 m in the vertical direction according to the thickness of the H₃b₁ sublayer. $314×340×40$ grids are finally determined.

The 3D training image is constructed using the object-based modeling method in commercial modeling software (Figure 11). The training image that can reflect the spreading characteristics of the sand body in the study area is established. The training image is a digital prototype model of the reservoir, which does not have to be faithful to the actual data but only an a priori geological concept. By comparing the training images (Figure 11) with the results simulated by the new algorithm (Figure 12) and the Snesim algorithm (Figure 13), it can be visually seen that the channel bars simulated by the new method are more geometrically complete and more similar to the training image built based on the a priori geological concept. The pattern of contact of sand bodies in space is consistent with the pattern of sand bodies’ architecture distribution in braided rivers.

According to the histogram of the statistical distribution of the thickness of the channel bar of the well data (Figure 14), it can be seen that the thickness of the channel bar in the work zone is mainly distributed in the range of 2 m-8 m, with a normal distribution trend. The statistical results of the simulation of the new method accounted for about 67.8% of 2-8 m, 0-2 m, 8-12 m, and more than 12 m accounted for about 10%, respectively. The trend of the normal distribution of the program has a high agreement with the statistical results of the well data. And the statistical results of Snesim method show a decreasing trend in general. The peak thickness is 0-2 m, which does not match with the actual distribution. It can be seen from Table 2 that the sample variance of the model simulated with the Snesim method is significantly smaller than the sample variance of the well data. And the sample variance of the new method is close to the sample variance of the well data.

From the table of variogram fitting parameters (Table 3), it can be seen that the new method is closer to the training image simulation results in terms of nugget, major, and minor range compared to the Snesim method simulation results. It can be visualized from the variogram model (Figure 15) that the training image fits better with the new method. The new method can better characterize the scale and morphology of the channel bars and characterize its integrity by better expressing the variation pattern of the length of the channel bars in the primary direction and the width of the channel bars in the minor direction.

Figure 16 shows the histogram of the distribution of the volume of a single channel bar. The simulation results of the new method show that the distribution probability values show a trend of increasing and then slowly decreasing as the volume of the channel bar increases. The distribution is mainly in the range of 50-550. Most of the volume of the channel bars simulated by the Snesim method are less than 300 grids. The number of simulated channel bar is large, but the size of the channel bar is too small.

The above study shows that both geostatistical modeling methods have some graphical reproduction capability. A careful analysis of each simulation implementation reveals that the new algorithm has a significant advantage over Snesim in characterizing the geometry, size, and stacking pattern of the channel bar. Snesim is also able to reproduce the training image information to some extent, but the simulation results are more discrete and do not retain the geometry of the channel bar well.

Braided river reservoir is a very important oil and gas reservoir type with complex internal inhomogeneity. To effectively characterize its internal inhomogeneity, fine characterization of architectural elements such as channel bars and braided channels is required. Traditional modeling methods have certain shortcomings in characterizing the internal architecture of the braided river reservoir. For example, it is difficult to characterize the scale of channel bars and to reproduce the internal architecture patterns of braided rivers by means of a sequential indication simulation method based on the variogram. The object-based modeling method can better portray the size of the channel bars, which can provide a strong foundation for the subsequent simulation of the silt layer. However, it is difficult to portray the internal construction pattern of braided river reservoir and the stacking patterns of composite channel bars. Snesim is a multipoint geostatistical modeling method. It can better satisfy the condition data and match the internal configuration pattern of braided river to some extent. However, the size of the channel bars is difficult to portray reasonably. This paper proposes a new modeling method for braided river reservoir architecture, which better solves the problems of geometry, scale, and conditionality of channel bars. Since the simulated channel bar is a complete object, it provides a good constraint for further simulation of the drop silt layer. Currently, the simulation is mainly for few wells, i.e., one well is drilled to encounter only one channel bar. If multiple wells are drilled to encounter a single channel bar, more constraints can be added. If the well spacing is larger than the plane size of the single channel bar, it can be designed as a composite of multiple channel bar. If the well spacing is smaller than the plane size of a single channel bar, the top and bottom surfaces of the channel bar are modified with the help of a face-based modeling strategy to achieve conditionality.

Due to the nature of this research, participants of this study did not agree for their data to be shared publicly, so supporting data is not available.

On behalf of all authors, the corresponding author states that there is no conflict of interest.

This research work was funded by the National Science and Technology Major Special Projects (2016ZX05027-004-003). In particular, we would like to thank Prof. Shaohua Li from Yangtze University for his guidance and suggestions.